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ABSTRACT 


The  free  electron  laser  (PEL)  is  proposed  to  meet  the  Navy’s  need  for  a  speed- 
of-light  high  energy  laser  weapon  capable  of  engaging  a  variety  of  targets  including 
anti-ship  cruise  missiles,  small  boats,  and  theater  ballistic  missiles.  A  key  attribute 
of  FELs  is  good  optical  beam  quality;  in  other  words,  they  operate  in  only  a  few  of 
the  lowest-order  transverse  Gaussian  modes.  For  weapons  applications,  a  good  mode 
quality  is  desired  because  it  delivers  the  highest  intensity  on  target  ensuring  a  high 
level  of  lethality.  A  few  higher-order  modes  can  arise  from  the  interaction  of  the 
electron  beam  with  the  optical  beam,  or  from  misalignments  of  the  electron  beam  or 
resonator  mirrors.  High  intensity  on  EEL  optics  can  lead  to  mirror  distortion  due  to 
heating  and  insufficient  cooling  of  the  mirror  substrate.  Mirror  distortions,  including 
astigmatism,  can  cause  higher-order  modes  to  appear  affecting  EEL  performance. 
Therefore,  it  is  important  to  quantify  these  higher-order  modes  because  doing  so 
uniquely  identifies  the  optical  field  and  may  allow  for  corrective  optics  to  single  out 
the  best  modes  for  EEL  lethality. 

This  thesis  will  review  free  electron  laser  theory,  and  for  the  first  time  develop 
analytical  solutions  to  quantify  Hermite-Gaussian  higher-order  modes,  develop  a  di¬ 
agnostic  for  modal  analysis,  and  determine  the  tolerance  limits  on  mirror  distortions. 
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I. 


MOTIVATION 


Previous  Naval  Postgraduate  School  (NPS)  classroom  studies  have  shown  that 
the  PHALANX  close-in  weapon  system,  used  onboard  US  Navy  ships  as  a  last  tier 
self-defense  weapon,  typically  destroys  an  enemy  missile  at  a  distance  of  a  few  hundred 
meters  from  the  ship.  At  that  distance,  NPS  classroom  studies  have  also  shown  a 
large  probability  that  debris  will  cause  collateral  damage  to  the  ship.  A  key  factor 
in  determining  debris  field  collateral  damage  is  time.  An  anti-ship  cruise  missile 
(ASCM)  typically  has  a  velocity  of  1200  m/s.  If  it  is  detected  at  a  nominal  horizontal 
distance^  of  15  km,  then  the  ship  has  approximately  13  seconds  to  track  and  engage 
the  missile.  The  probability  of  collateral  damage  stems  from  the  fact  that  it  takes 
several  PHALANX  rounds  to  destroy  the  missile.  So,  Naval  interest  in  high  energy 
lasers  is  clear:  by  engaging  a  target  with  a  speed-of-light  high  energy  laser  weapon,  the 
distance  and  time-to-engage  increases  considerably  and  the  probability  of  collateral 
damage  to  the  ship  is  reduced. 

In  order  to  “kill”  a  missile,  about  one  liter  of  missile  material  must  be  de¬ 
stroyed.  This  will  cause  sufficient  structural  damage  resulting  in  missile  breakup. 
The  amount  of  energy  required  for  material  ablation  is  on  the  order  of  one  megawatt 
over  a  surface  area  of  100  cm^  for  a  few  seconds.  If  an  ASCM,  like  the  one  described 
above,  is  detected  at  15  km  and  engaged  by  a  high  energy  laser,  it  would  be  destroyed 
in  roughly  five  seconds  at  a  distance  of  9  km  from  the  ship  reducing  the  probability 
of  collateral  damage  considerably  giving  the  ship  a  large  safety  margin  in  time  and 
distance. 

Several  U.S.  Navy  laser  programs  are  currently  investigating  the  feasibility  of 
employing  high  energy  lasers  on  board  ships.  The  free  electron  laser  (PEL)  possesses 

^Horizontal  distance  as  a  function  of  height  of  eye  is  given  by  (nautical  miles)  = 
1.169v^Height  of  eye  (ft) 
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many  attributes  that  are  in-line  with  the  the  Navy’s  vision  of  an  all-electric  ship  with 
no  munitions: 

•  High  average  power 

•  Can  be  designed  over  a  range  of  wavelengths  (IR  through  X-rays) 

•  Good  optical  beam  quality 

•  High  operational  reliability 

•  No  hazardous  chemicals  required  for  lasing  medium 

•  Virtually  limitless  magazine  depth 

•  Capable  of  engaging  a  variety  of  targets  including  ASCMs,  small  boats,  and 
theater  ballistic  missiles 

A  key  attribute  of  FELs  is  good  optical  beam  quality;  in  other  words,  they  operate 
in  only  a  few  of  the  lowest-order  Gaussian  modes.  A  key  attribute  of  FELs  is  good 
optical  beam  quality;  in  other  words,  they  operate  in  only  a  few  of  the  lowest-order 
transverse  Gaussian  modes.  For  weapons  applications,  a  good  mode  quality  is  desired 
because  it  delivers  the  highest  intensity  on  target  ensuring  a  high  level  of  lethality.  A 
few  higher-order  modes  can  arise  from  the  interaction  of  the  electron  beam  with  the 
optical  beam,  or  from  misalignments  of  the  electron  beam  or  resonator  mirrors.  High 
intensity  on  EEL  optics  can  lead  to  mirror  distortion  due  to  heating  and  insufficient 
cooling  of  the  mirror  substrate.  Mirror  distortions,  including  astigmatism,  can  cause 
higher-order  modes  to  appear  affecting  EEL  performance.  Therefore,  it  is  important 
to  quantify  these  higher-order  modes  because  doing  so  uniquely  identifies  the  optical 
field  and  may  allow  for  corrective  optics  to  single  out  the  best  modes  for  EEL  lethality. 

This  thesis  will  review  free  electron  laser  theory,  and  for  the  first  time  develop 
analytical  solutions  to  quantify  Hermite-Gaussian  higher-order  modes,  develop  a  di¬ 
agnostic  for  modal  analysis,  and  determine  the  tolerance  limits  on  mirror  distortions. 
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II. 


FREE  ELECTRON  LASER  SYSTEM 


The  free  electron  laser  (PEL)  system  consists  of  many  different  components 
each  of  which  involves  very  interesting  physics.  Figure  1  shows  a  basic  systems  di¬ 
agram  of  the  operating  EEL  at  the  Thomas  Jefferson  National  Accelerator  Facility 
(JLab).  Let  us  walk  through  the  operation  of  JLab’s  EEL.  A  pulsed  beam  of  electrons 
in  a  vacuum  is  produced  using  a  photoinjector.  The  beam  of  electrons  are  accelerated 
to  near  the  speed  of  light  in  the  superconducting  accelerator  using  radio-frequency 
(RE)  energy.  Highly  relativistic  electrons  are  then  “wiggled”  in  the  undulator  pro¬ 
ducing  laser  light  where  it  is  amplified  by  extracting  energy  from  the  electron  beam 
within  the  optical  resonator.  The  electron  beam’s  remaining  energy  is  recycled  by 
having  them  enter  the  superconducting  accelerator  out  of  phase  with  the  RE  fields. 
Finally,  the  electron  beam  is  directed  to  a  beam  dump. 

At  the  heart  of  the  EEL  system  are  relativistic  “free”  electrons,  i.e.,  not  bound 
to  any  atom  or  molecule.  The  electrons  travel  in  a  vacuum  much  like  in  a  cathode  ray 
tube.  Conventional  lasers  operate  at  specific  frequencies  because  electrons  can  make 
transitions  only  between  discrete  energy  levels.  The  electrons  of  an  EEL  are  forced 
to  vibrate  through  a  spatial  periodic  magnetic  field.  The  frequency  of  vibration  can 
be  changed  by  adjusting  the  magnetic  field  or  by  changing  the  speed  of  the  electrons 
which  in  turn  changes  the  laser  light  frequency. 

In  1970  Dr.  John  Madey  of  Stanford  University  proposed  what  he  termed 
a  “free  electron  laser”  using  highly  relativistic  electrons  to  produce  radiation  in  a 
magnetic  undulator.  Six  years  later  Madey  succeeded  in  demonstrating  gain  with  a 
free  electron  laser  with  a  24  MeV  electron  beam  and  5  m  long  undulator  [2]. 

As  an  illustrative  example  of  an  operating  EEL,  Table  I  shows  the  current 
parameters  at  JLab.  At  the  time  of  this  writing,  JLab  has  the  most  powerful  free 
electron  laser  in  the  world  at  10  kW.  The  wavelength  is  in  the  infrared  around  2-4 
//m. 


3 


^  recycled  electrons 

t 

electron  beam 

optical 
cavity  | 

mirrors  1 

)  lastr  light 

* 

beam 

dump 

source  superconducting 

wiggler  for  conversion 

accelerator  essential  for 

of  electron  energy 

recycling  electron  energy 

into  light 

Figure  1.  Free  electron  laser  systems  diagram  highlighting  major  components  from 

[!]• 


A.  HOW  DOES  IT  WORK? 

For  an  FEL  to  work  we  need  [4]: 

•  a  beam  of  relativistic  electrons, 

•  co-propagating  in  an  optical  field, 

•  through  a  series  of  spatially  periodic  magnets. 


The  undulator  establishes  electron  oscillations  that  are  transverse  to  its  overall  longi¬ 
tudinal  motion.  In  a  reference  frame  moving  along  the  FEL  axis  with  the  electrons, 
it  would  appear  that  they  are  simply  oscillating  and  radiating  in  all  directions.  In  the 
laboratory  frame  the  radiation  is  confined  to  a  narrow  forward  cone  whose  angle  is  a 
function  of  the  electron’s  speed,  and  is  known  as  magnetic  Bremsstrahlung  radiation 
[2].  The  radiation  emitted  is  highly  Doppler  shifted  due  to  the  relativistic  nature  of 
the  electrons  and  the  wavelength  is  given  by 


A 


A, 


(K^  +  1) 

'  272 


(11.1) 
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Optical  Resonator 

Output  power 

10  kW 

Length 

32m 

Wavelength 

l.l//m 

Waist  radius 

0.6  mm 

Mode  radius  at  mirror 

1  cm 

Quality  factor 

25  (4%  outcoupling) 

Rayleigh  length 

0.5  m 

Max.  intensity  permissible  on  mirrors 

200  kW/cm^ 

Electron  Beam 

Energy 

111  MeV 

Bunch  charge 

135  pC 

Radius 

0.2  mm 

Peak  current 

340  A 

Pulse  length 

0.12  mm 

Undulator 

Period,  Ao 

8  cm 

Length 

2.3 

Number  of  periods 

29 

Undulator  parameter,  K 

0.55 

Table  I.  JLab  PEL  Parameters  from  [3]. 


where  A  is  the  wavelength  of  the  optical  radiation,  Ao  is  the  wavelength  of  the  undu- 
lator  magnetic  field,  7  is  the  Lorentz  factor,  K  is  a  dimensionless  parameter  that  is 
proportional  to  the  undulator  magnetic  field  strength.  We  will  derive  this  equation 
in  the  next  chapter. 

The  initial  optical  radiation  given  off  by  the  oscillating  electrons  is  known  as 
spontaneous  emission  in  the  vernacular  of  conventional  laser  physics.  The  radiation 
from  spontaneous  emission  is  incoherent,  meaning  that,  the  wave  packets  emanating 
from  each  radiating  electron  are  not  in  phase  with  adjacent  wavepackets  since  each 
electron  is  randomly  dispersed  within  the  electron  beam.  Coherence  evolves  from  this 
cacophony  through  mode  competition  with  many  passes  through  the  undulator.  The 
process  of  mode  competition  will  be  covered  later. 
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At  resonance,  the  electric  field  felt  by  the  electrons  oscillates  at  the  same 
frequency  and  in  the  same  plane  as  the  electrons  oscillating  through  the  periodic 
undulator  magnetic  field.  If  the  electric  field  and  the  electron’s  transverse  motion  is 
parallel,  electrons  will  lose  energy  and  slow  down  since  work  is  being  done  by  them. 
If  the  electric  field  and  the  electron’s  transverse  motion  is  anti-parallel,  electrons  will 
gain  energy  and  speed  up  since  work  is  being  done  on  them.  The  energy  they  lose  or 
gain  comes  from  the  optical  radiation  field.  The  effect  of  electrons  losing  and  gaining 
energy  within  each  optical  wavelength  is  a  longitudinal  bunching  of  electrons  as  some 
slow  down  and  others  speed  up  in  each  wavelength.  The  bunched  electrons  then  emit 
radiation  that  is  in  phase  and  coherent.  The  intensity  of  the  coherent  radiation  is 
much  more  powerful  than  the  incoherent  radiation  [2]. 

Two  key  performance  metrics  of  the  FEL  can  now  be  introduced:  extrac¬ 
tion  and  gain.  The  fractional  increase  in  optical  radiation  power  in  one  pass  of  the 
undulator  is  defined  as  gain, 

P  —  P 

G  =  (II.2) 

where  the  optical  power,  P,  is  measured  at  the  end  of  the  undulator  and  Pg  is  the 
power  at  the  beginning  of  the  undulator.  As  the  FEL  converts  electron  kinetic  energy 
to  optical  radiation  energy,  the  fractional  energy  extracted  from  the  electron  beam  in 
one  pass  of  the  undulator  is  defined  as  extraction,  77, 


Optical  power 

Tj  =  — - 

Initial  electron  beam  power 


(11.3) 


The  optical  gain  must  be  sufficient  to  overcome  loss  and  outcoupling.  A  larger  ex¬ 
traction  means  the  FEL  requires  less  average  current  to  achieve  high  power.  Typical 
values  for  extraction  and  gain  are  2%  and  30%,  respectively. 


B.  CONFIGURATIONS 

When  the  undulator  is  bounded  within  two  reflective  mirrors,  the  configura¬ 
tion  is  known  as  an  optical  oscillator  or  resonator.  The  optical  radiation  field  that 
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couples  to  the  electron’s  motion  can  manifest  itself  through  spontaneous  emission. 
This  radiation  is  collected  within  mirrors  as  shown  in  Figure  2.  A  partially  reflective 
outcoupling  mirror  provides  the  optical  feedback  back  into  the  undulator.  Coherence 
develops  through  many  passes  of  the  optical  held.  The  source  of  the  optical  radi- 


Figure  2.  Free  electron  laser  system  in  an  oscillator  conflguration  from  [5]. 


ation  held  can  also  be  an  external  laser  in  the  case  of  an  amplifler.  An  amplifler 
conflguration  makes  use  of  a  seed  laser  and  typically  a  longer  undulator.  The  longer 
undulator  provides  more  gain  and  extraction  in  a  single  pass  since  only  one  pass  is 
used.  Coherence  is  established  by  the  seed  laser.  Figure  3  shows  a  typical  amplifler 
conflguration. 


Electron  Beam 


B  field 


Accelerator  ■  fc--"  ■■■titttitiwt- 


Laser  Light 


^See^Lase^^^ 


Figure  3.  Free  electron  laser  system  in  an  amplifler  conflguration  from  [5]. 
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C.  COMPONENTS 
1.  Source 


Photoinjector  schematic 

RF 

m 


.'  ■ill  I 
'■rillatOf 


Photocathode 


r- 

hijii-.tor  RF  Linac 

i 

■ 
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0.5- 1  MeV 


7-10  MeV 


Figure  4.  Schematic  of  a  typical  photoinjector  from  [6]. 


An  FEL  requires  high  peak  current  in  order  for  coherent  radiation  to  grow  in 
amplitude.  The  source  for  free  electrons  is  a  photoinjector.  A  drive  laser,  usually 
a  solid-state  laser,  is  pulsed  at  a  photocathode  target.  Through  the  photoelectric 
effect,  electrons  are  literally  blasted  from  the  surface  of  a  cathode  and  accelerated  in 
an  electric  field  to  the  linear  accelerator.  Since  the  drive  laser  fires  at  regular  intervals 
(on  the  order  of  nanoseconds),  electrons  travel  in  pulses.  Electrons  at  the  terminus 
of  the  photoinjector  have  energies  from  5  to  10  million  electron-volts  (MeV). 

2.  Accelerator 

The  linear  accelerator  (see  Eigure  5)  has  two  functions:  (1)  accelerate  electrons 
from  the  photoinjector  to  its  final  velocity  (energy  100  MeV),  and  (2)  recover 
energy  from  recycled  electrons.  The  electrons  are  accelerated  as  they  travel  through 
the  cavities  in  phase  with  the  radio- frequency  (RE)  fields.  At  JLab,  electron  bundles 
from  the  photoinjector  are  accelerated  using  a  superconducting  RE  linear  accelerator 
to  minimize  waste  heat.  Energy  is  recovered  from  recycled  electrons  and  given  back 
to  the  RE  field  since  they  reenter  the  accelerator  out  of  phase  with  the  RE  field. 
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Figure  5.  JLab  linear  accelerator  from  [7] 


3.  Undulator 

The  undulator  is  a  collection  of  spatially  alternating  permanent  magnets  that 
causes  the  electrons  to  wiggle  transverse  to  their  longitudinal  motion.  The  transverse 
motion  then  couples  to  the  optical  electric  field  enabling  work  to  be  done  on  or  by  the 
optical  radiation  field.  The  magnets  can  be  orientated  linearly  or  helically  as  shown 
in  Figure  6.  Linear  undulators  produce  linearly  polarized  light  whose  electric  field 
oscillates  in  the  same  plane  as  the  oscillating  electrons.  Helical  undulators  produce 
circularly  polarized  light  whose  electric  field  rotates  along  with  the  electrons  enabling 
a  coupling  to  occur. 

4.  Beam  Dump 

After  energy  is  recovered  from  the  recycled  electrons  in  the  linear  accelerator, 
the  electron  beam  is  steered  into  a  beam  dump  (see  Figure  7).  A  beam  dump  is 
nothing  more  than  a  piece  of  notched  metal  that  will  absorb  the  remaining  kinetic 
energy  of  the  electron  beam  (on  the  order  of  several  MeV).  At  at  an  average  electron 
beam  current  of  one  ampere,  about  one  megawatt  must  be  absorbed.  It  is  water- 
cooled  on  the  non-vacuum  side. 
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Figure  6.  (a)  Helical  undulator  resulting  in  circularly  polarized  light,  (b)  Linear 
undulator  resulting  in  linearly  polarized  light  from  [8]. 

5.  Optics 

As  mentioned  in  Table  I,  the  maximum  permissible  intensity  on  the  resonator 
mirrors  is  on  the  order  of  200  kW/cm^  to  prevent  mirror  damage  to  the  mirror’s 
dielectric  layers.  In  adequate  or  nonuniform  mirror  cooling  may  also  warp  the  sapphire 
substrate  thereby  affecting  FEL  performance.  FEL  designs  with  a  short  Rayleigth 
length^  have  modes  that  expand  so  that  the  intensity  on  the  mirrors  is  lessened.  The 
entire  optical  system  includes  all  the  necessary  beam  directors  to  place  laser  energy 
on  a  target  miles  away  from  the  ship. 


^Rayleigh  length  is  the  distance  the  mode  has  to  travel  longitudinally  from  the  waist  so  that  the 
mode  area  doubles.  Rayleigh  length  will  be  discussed  in  greater  detail  in  Chapter  IV. 
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Figure  7.  Beam  dump  from  [9] 
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III. 


FREE  ELECTRON  LASER  THEORY 


We  will  now  develop  in  detail  the  mathematics  that  describes  the  operation 
of  the  free  electron  laser.  In  order  to  produce  laser  light  we  must  cause  electrons  to 
oscillate  through  a  Lorentz  force  interaction,  bunch  the  electrons  spatially  for  coher¬ 
ence,  and  finally,  have  sufficient  electron-to-optical  mode  coupling  for  amplification 
and  gain. 


A.  EQUATIONS  OF  MOTION 
1.  Relativistic  Lorentz  Force 

In  compact  covariant  notation,  the  relativistic  Lorentz  Force  law^  is 

e 


dp"  df/" 
=  m- 


=  — 
c 


dt'  dt'  c' 

where  p"  =  {jmc,  p)  is  the  four-momentum,  7  =  1/  jc?  is  the  Lorentz  factor, 

V  is  the  electron’s  velocity,  m  is  the  mass  of  the  electron,  =  (7c,  —7V)  is  the 
electron’s  four- velocity,  t'  =  t/'j  is  proper  time,  e  is  the  magnitude  of  the  electron’s 
charge,  c  is  the  speed  of  light,  and  is  the  field-strength  tensor  with  the  components 
of  the  electric  field,  E,  and  magnetic  field,  B.  Specifically, 

/ 


_ 


V 


0 

~Ex 

Ey 

-E, 

Ex 

0 

-B, 

By 

Ey 

B, 

0 

-Bx 

E, 

—  By 

Bx 

0 

(III.2) 


y 


Let  us  describe  the  electron’s  trajectory  in  a  helical  undulator  whose  magnetic  field 
along  2;  is 

Bu  =  B  o{coskoZ,sinkoZ,0)  (III. 3) 

where  Bo  is  the  magnetic  field  amplitude,  ko  =  2n/Xo  is  the  undulator  wavenumber, 
and  Ao  is  the  undulator  wavelength,  i.e.,  the  distance  for  one  complete  cycle  of  the 


^[10],  Chapter  11 
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undulator  magnetic  field.  The  corresponding  electric  and  magnetic  field  of  the  optical 
radiation  field  is: 

Er  =  £'o(cos'0,  —  sin?/),  0) 

Br  =  Eo{sm  ?/),  cos  ?/),  0)  (III. 5) 

where  'll)  =  kz  —  ut  +  (f),  k  is  the  optical  wavenumber,  u  =  kc  is  the  optical  angular 
frequency.  Eg  is  the  field  amplitude  in  cgs  units,  and  (f)  is  the  optical  phase. 

If  we  combine  Equations  (HI. 3),  (HI. 4),  and  (HI. 5),  Equation  (HI. 2)  becomes 

^  0  —Eg  cos  (f)  Eg  sin  cf)  0 

Eg  COS  ?/)  0  0  Bo  sin  kgZ  +  Eg  cos  ^ 

—Bo  sin?/)  0  0  —  Bo  cos  A:o.2:  —  Bo  sin 

0  —Bo  sin  kgZ  —  Eg  sin  xp  Bg  cos  kgZ  +  Eg  sin  xp  0 


pa/3  _ 


V 


Using  Equations  (III.l),  (III-6),  and 


/3  =  - 

c 

7  = 


(111.6) 

(111.7) 

(111.8) 


where  (3  is  dimensionless  velocity  and  7  is  the  Lorentz  factor.  We  can  write  down  the 
spatial  components  of  the  force  from  Equation  (III.l), 

=  -  —  [Eg(l  -  l5z){cosxp,-smxp,Q)  +  l5zBg{-siTikgZ,coskgZ,Q)] 

at 


d(7/5. 


me 

e 


g 

= - [EgifixCOSxp  -  l^ySmxp)  ^  BgifixSiwkgZ  -  P^yCOSkgZ)]  (III. 9) 

at  me 

where  we  have  made  use  of  f3±  =  I3xi  +  I3yj.  Since  /5^  1  ^  (1  —  Pz)  — >■  0,  we  then  drop 

the  terms  proportional  to  Bo(l  —  /3z)  compared  to  /5^Bo.  Integrating  the  transverse 
equation  gives  (assuming  perfect  injection^) 

— gBq  — cBqXq 

/3  I  = - ^(cos  kgZ,  sin  kgZ,  0)  = - ^(cos  kgZ,  sin  kgZ,  0)  (III. 10) 

jme^kg  zTT'jme^ 


^Perfection  injection  means  that  the  electrons  enter  the  undulator  at  precisely  the  correct  angle 
such  that  its  motion  is  exactly  given  by  Equation  (III.  10)  and  constants  of  integration  are  zero. 


\ 


J 
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We  can  clean  this  equation  up  a  little  if  we  identify  the  dimensionless  undulator 
parameter,  K  =  eBo\o/{‘2-nrnc^)  to  write: 

-K 

I3j_  =  - (cos  koZ,  sin  koZ^O)  (III.ll) 

7 

We  can  now  derive  the  energy  exchange  equation  between  the  electrons  and 
the  optical  field.  With  Equation  (III. 2),  dt'  =  dt/q  and  =  Ejc  =  'jmc,  from  rela¬ 
tivistic  kinematics^,  we  can  write  down  the  zeroth  component,  or  energy  component, 
of  Equation  (III.l).  Eor  a  =  0,  we  have 

dr  c 


so  that 


M/  , 

'y—i'ymc 

'dV' 


--(0, -E)  •  (yc,  -yv) 
c 


and 


d  ,  , 

= - E-v 

'dV'  '  c 

giving 

^  /3-E  (III.12) 

dt  me 

Equation  (III.  12)  describes  how  the  the  electron  exchanges  energy  with  the  optical 

field.  The  electron’s  total  energy,  £  =  is  a  function  of  its  rest  mass  and  its 

speed  (through  7).  The  undulator  sets  up  electron  motion  in  the  transverse  directions 
otherwise  /3  ■  E  would  be  zero. 

Substituting  Equations  (HI. 4)  and  (III.ll)  into  (III. 12),  we  arrive  at 

eKE 

7=  cos(koZ  +  ij)  (III.13) 

jmc 


^see  [11],  Chapter  1 
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Recall  that  'll)  =  kz  —  iot  +  (j).  If  we  define 

(  =  (ko  +  k)z  —  ujt  (111.14) 

as  the  electron  phase,  we  have 

eKE 

7=  cos ((  +  (/>)  (III.15) 

jmc 

We  have  arrived  at  an  expression  that  describes  energy  exchange  between  the  electron 
and  optical  field.  Electron  phase,  (,  is  another  way  of  describing  the  electron’s  position 
along  the  undulator  axis  with  respect  to  the  co-propagating  optical  wave.  Notice  that 
for  cos(C  +  (j))  >  0,  the  electron  energy  increases  (absorption);  for  values  less  than 
zero,  the  electron  loses  energy  (stimulated  emission). 

2.  Resonance  Condition  and  Dimensionless  Variables 
Introduced 

Since  the  electron  phase,  is  proportional  to  the  electron  position,  2;,  let  us 
define  the  dimensionless  electron  phase  velocity  as 

v=CL/c  (III.16) 

where  L  is  the  length  of  the  undulator.  We  can  normalize  the  evolution  of  the 
electron’s  motion  to  the  length  of  the  undulator  by  defining  a  dimensionless  time, 
T  =  ct I L.  It  is  convenient  when  possible  to  normalize  our  equations  on  scales  that 
are  relevant  to  the  physics  being  explored.  For  example,  longitudinal  length  can  be 
normalized  to  the  length  of  undulator.  Thus,  r  takes  on  values  from  0  — >■  1  from  the 
beginning  of  the  undulator  to  the  end  of  the  undulator.  In  the  following  sections,  we 
will  normalize  transverse  directions,  optical  field  amplitude,  and  current  density. 
Returning  to  phase  velocity,  v  can  also  be  expressed  as 

D=C  (III.17) 

where  (.?.)  =  d{. .  .)/dT  is  a  derivative  with  respect  to  dimensionless  time  r.  Phase 
velocity,  v,  is  a  measure  of  how  fast  the  electron  is  moving  longitudinally  within  an 
undulator  period.  We  can  now  derive  the  resonance  condition. 
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The  resonance  condition  results  in  optimum  energy  exchange  between  the  elec¬ 
tron  and  the  optical  mode.  This  occurs  when  one  wavelength  of  light  passes  over  an 
electron  as  the  electron  travels  through  one  undulator  period.  PEL  resonance  occurs 
when  =  0.  Taking  the  time  derivative  of  the  electron  phase  gives  us 

C  =  (ko  +  k)z  —  ujt 

C  =  {ko  +  k)PzC  —  uj  (III. 18) 


Substituting  this  expression  into  the  definition  of  electron  phase  velocity  gives  us 


^  =  [(ko  +  k)l3zC  —  uj]  — 

c 


(III.19) 


For  =  0,  we  have 


k  = 


koPz 


so  that  A  =  Ao 


(1  -  A) 

(1  -  A) 


A 


(III.  20) 


where  we  have  used  k  =  u/c.  So  now  let  us  find  an  expression  for  j5z  in  terms  of 
physical  parameters  of  the  EEL.  From  Equation  (III.ll),  /5j_  =  K/y,  and  the  definition 
of  the  Lorentz  factor,  we  have 

1 

7  = 


7^  = 


7^  = 


\/r^ 

1 

1  -  /5i  -  /5f 


1  _  ^  _  /52 
^2  Pz 


(III.21) 


Solving  for  the  electron  velocity  gives  us 


/5,  =  Wl--(K2  +  i) 


T 


(III.22) 
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For  highly  relativistic  electrons  employed  in  an  FEL,  7  is  typically  100,  so  we  can 
make  use  of  the  binomial  expansion"^  giving  us 

A  =  + 

/5.  ~  1-^(K2  +  1)  (IIL23) 

Returning  to  Equation  (III 20),  we  then  have. 


A 


A, 


+  1) 

’  272 


(IIL24) 


Here  is  an  expression  for  the  wavelength  of  the  optical  mode  based  on  physical  pa¬ 
rameters  that  we  can  change  (Ao,  -R,  7):  a  clear  example  of  the  EEL’s  tunability.  It 
is  now  understandable  why  the  wavelength  of  the  light  is  so  much  smaller  (by  a  factor 
of  1/7^)  than  the  undulator  period.  This  is  physically  plausible  since  the  relativistic 
nature  of  the  electrons  causes  them  to  see  the  undulator  Lorentz-contracted  by  a 
factor  of  7  and  Doppler-shifted  light  from  the  electrons’  frame  is  then  transformed  to 
the  lab  frame  reducing  the  wavelength  by  another  factor  of  7.  Electrons  see  resonant 
forces  from  undulator  and  the  optical  field.  The  end  result  is  the  optical  field  is  at  a 
much  higher  frequency  then  the  undulator’s  spatial  frequency. 


3.  Pendulum  Equation 

We  are  now  in  a  position  to  recast  the  equations  in  the  previous  sections  into 
a  form  known  as  the  pendulum  equation.  In  this  pendulum  description,  phase-space 
plotting  becomes  a  natural  and  highly  illustrative  way  of  describing  the  electrons’ 
motion  and  interaction  with  the  optical  mode.  Starting  with  the  phase  velocity,  we 


"^The  binomial  expansion  for  small  x  is  given  by 


{l  +  x)'^  =  l  +  nx+  . 

see  [12]  pp.  14  for  more  details 
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have 


(111.25) 

(111.26) 


p  =  [{ko  +  k)(5z-k]L=C, 

p  =  (kg  k)  f3z  L 

So  let  us  now  calculate  an  expression  for  /5^.  From  Equation  (III.  23), 

A  =  +  (111.27) 

/5.  =  -AK‘^  +  1)  (III.28) 

But  near  resonance  we  have  so  Equation  (III. 28)  becomes 

„  _  7(A-"  +  l) 

Pz  9 

7  7^ 

A  2  A 

/5.  =  -V  (III.29) 

7  ^0 

Returning  to  Equation  (III. 26),  substituting  Equation  (III. 29),  and  using  k  =  27r/A, 
we  arrive  at 

p  = 

p  = 

We  can  make  the  following  substitutions  and  approximations:  since  the  undulator 
period  is  much  larger  than  the  optical  wavelength  (i.e,  Xg  »  A),  and  since  the  first 
term  is  smaller  than  the  second  by  a  factor  X/Xg,  the  first  term  is  negligible.  We  define 
|a|  =  AnNeK\E\L/AmA  as  the  dimensionless  optical  field  amplitude,  L  =  NXg,  the 
length  of  the  undulator,  dt  =  Ldr /c,  and  substituting  Equation  (III. 15),  we  finally 
arrive  at 

O  °° 

|a|  cos(C  +  ^)  =C  (III. 32) 

The  pendulum  equation  describes  the  motion  of  the  electrons  in  phase  space,  (Cw)- 


47rA  7  47r  7 


(111.30) 

(111.31) 
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4.  Phase  Space 

Phase  space  diagrams  are  used  in  many  physical  problems.  For  example,  in  a 
simple  pendulum,  the  angular  displacement,  9{t),  and  its  angular  velocity,  z^(t),  can 
be  thought  of  coordinates  of  a  two-dimensional  phase  space.  As  the  motion  of  the 
pendulum  evolves,  the  point  given  by  (6,  v)  will  move  along  a  path  in  phase  space. 
For  different  initial  conditions,  the  motion  will  be  described  by  different  paths.  Phase 
space  plots  can  be  divided  into  two  categories:  closed  orbit  or  open  orbit  paths. 

For  the  pendulum,  closed  orbit  paths  result  when  the  pendulum  is  bound 
about  the  stable  fixed  point.  Open  orbit  paths  occur  when  there  is  a  large  enough 
angular  velocity  to  cause  the  pendulum  to  go  over  the  top. 


Pendulum  Phase  Space  Paths 


0 


Figure  8.  Various  phase-space  paths  for  a  simple  pendulum.  The  thin  green  curves 
are  examples  of  closed  orbits.  The  dotted  red  curves  are  examples  of  open  orbits. 
The  separatrix  is  the  heavy  dark  curve. 


Figure  8  illustrates  the  phase  space  paths  of  a  simple  pendulum  for  10  different 
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initial  conditions.  The  phase  space  paths  are  given  by® 

=  vl  +  2(cos  9  —  cos  9o)  (111.33) 

The  separatrix  is  found  when  9o  =  tt  and  =  0  in  Equation  (III. 33)  giving  us 

=  2(cos0  + 1)  (III.34) 

The  separatrix  separates  locally  bounded  motion  from  locally  unbounded  motion  and 
passes  through  the  unstable  fixed  point  {6o  =  tt,  =  0). 

Equation  (III. 32)  is  equivalent  to  the  simple  pendulum  if  |a|  and  (f)  are  constant, 
as  is  approximately  the  case  for  low  gain  in  the  EEL.  The  electron  phase,  can  be 
viewed  as  the  electron’s  position  within  an  undulator  period  while  a  change  in  phase 
velocity  is  proportional  to  the  electron’s  energy  (recall  z^ocT/q).  The  dimensionless 
optical  field  amplitude,  |a|,  is  a  measure  of  the  field  strength.  A  larger  |a|  causes 
the  phase  space  paths  to  evolve  much  faster  (analogous  to  a  larger  ^/gjl  for  the  real 
pendulum).  A  larger  |a|  also  causes  the  separatrix  to  increase  in  size  (along  p).  The 
EEL  separatrix  is  given  by  [13]: 

=  2|a|[l  +  sin(((  +  ^)]  (111.35) 

Eigure  9  shows  the  evolution  of  1000  electrons  in  phase  space  all  with  an  initial 
=  0  (on  resonance)  and  |a|  =  1.  Notice  that  half  the  electrons  gained  energy  from 
the  optical  field  (increasing  values  of  v)  and  half  the  electrons  lost  energy  to  the 
optical  field  (decreasing  values  of  p).  Notice  how  the  electrons  are  becoming  bunched 
along  (.  Since  there  are  an  equal  number  of  electrons  gaining  and  losing  energy,  the 
overall  gain  is  zero.  Gain  is  the  fractional  change  in  the  power  of  the  optical  field 

^Dimensionless  pendulum  equation  given  by  8  =  —  sin  6.  Multiply  both  sides  by  v  to  get  v  v= 
—v&mO,  which  is  equivalent  to  d[u‘^/2  —  cosd]/dt  =  0,  giving  us  Equation  (III. 33). 


21 


0.1 


0 

0.1 

0 

-7c/2  t;  3jc/2  0  t  1 

Figure  9.  Phase  space  evolution  for  1000  electrons  with  all  with  initial  i/  =  0. 


Gain 
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over  a  pass  through  the  undulator.  In  symbols,  we  write: 


G 

and 

P 


P-Pq 

Po 


(III.36) 


(III.37) 


where  the  power,  P  is  measured  at  the  end  of  the  undulator  and  Pq  is  the  power  at 
the  beginning  of  the  undulator.  We  will  discuss  gain  in  detail  in  the  next  section. 

In  Figure  10,  we  start  out  with  an  initial  Vq  =  1-5  and  |a|  =  5.0.  There  is 
an  imbalance  in  the  number  of  electrons  that  are  gaining  and  losing  energy  from  the 
optical  field.  In  this  particular  case,  there  is  a  clear  bunching  of  electrons  near  C  =  ^ 
resulting  in  gain. 


B.  ELECTRON  BEAM  AND  OPTICAL  FIELD  INTER¬ 
ACTION 

In  the  previous  section,  we  derived  the  FEL  pendulum  equation  which  tells  us 
the  phase  space  evolution  of  electrons  in  an  undulator  period.  Phase  space  evolutions 
that  ended  with  good  bunching  at  the  appropriate  phase  resulted  in  gain  (see  Figure 
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Figure  10.  Phase  space  evolution  for  1000  electrons  all  with  initial  =  1.5. 


10).  The  upper  picture  of  Figure  11  illustrates  bunching  in  an  optical  wavelength  for 


Figure  11.  In  the  upper  picture,  electrons  have  an  initial  =  0  which  results  in 
bunching  but  no  gain.  In  the  lower  picture,  electrons  have  a  an  initial  I'o  >  0  resulting 
in  good  bunching  and  positive  gain  from  [4]. 


=  0  or  on-resonance  case.  At  time  ti  the  electrons  start  out  evenly  spaced  within 
a  section  of  the  electron  beam  one  wavelength  of  light  long.  At  time  t2,  the  electrons’ 
positions  have  evolved  so  that  the  number  of  electrons  which  give  energy  to  the  optical 
field  (gain)  equals  the  number  of  electrons  which  take  energy  from  the  optical  field 
(absorption).  The  result  is  zero  net  gain  but  good  bunching.  In  the  lower  picture,  at 


23 


time  ti  evenly  spaced  electrons  have  an  initial  non-zero  phase  velocity,  i.e.  u  >  0.  At 
time  t2  the  electrons’  positions  have  evolved  so  that  there  are  more  electrons  giving 
energy  to  the  optical  field  than  absorbing  it.  The  result  is  positive  gain,  i.e.,  the 
optical  field  grows. 

1.  Maxwell’s  Equations 

Let  us  quantitatively  describe  the  interaction  between  the  electron  beam  and 
the  optical  field.  We  begin  with  Maxwell’s  Equations  in  cgs  units®: 


V  E 

=  47rp 

(III.38) 

VB 

=  0 

1  dB 

(III.39) 

V  X  E 

c  dt 

47r  IdE 

(III.40) 

V  X  B 

=  — J  + 

c  c  ot 

(III.41) 

where  p  is  the  charge  density,  E  is  the  electric  field,  B  is  the  magnetic  field,  and  J 
is  the  current  density.  Equation  (III. 39)  automatically  implies  that  B  is  the  curl  of 
some  other  vector.  A,  the  magnetic  vector  potential^.  We  then  have 


VB=O^B=VxA 


(III.42) 


If  we  substitute  Equation  (III. 42)  into  Equation  (III. 40)  we  get 

1  d 

V  X  E  =  —  —  (V  X  A) 
c  at 

1  d 

V  X  E  +  -  — (V  X  A)  =  0 
c  at 

Recall  from  elementary  vector  calculus,  if  the  curl  of  some  vector  is  zero,  then  that 
vector  is  equal  to  minus  the  gradiant  of  a  scalar®,  $,  the  electric  potential.  We  can 


®see  [10],  Chapters  6  and  11 

’^[14],  Chapter  5 

Tbid 
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then  write 


_  19A\  1 9A 

V  X  (  E  +  =  0  E+  =  -V$ 

c  at  j  c  ot 


Therefore, 

1  .^A 

E  =  -V$ - ^  (111.43) 

c  dt  ^  ’ 

Equations  (111.42)  and  (III. 43)  describe  the  electromagnetic  radiation  in  terms  of  their 

potentials,  A  and  $.  If  we  substitute  them  into  Equation  (III. 41)  we  arrive  at^ 


VxVxA  =  — J  +  -|-(  -V$  -  1 

c  cot 


V(V-A)- V'A 

a  1 

V  A  - 


47r 


-J  -  -V  — , 
:  c  ot 

A 

c  ot 


c  dt  J 
1  d^A 
9  dt‘^ 

47r , 


(III.44) 


Eor  completeness,  if  we  substitute  Equation  (III. 43)  into  Equation  (III. 38)  we  get 


-v$- 

c  dt 
1  d 


=  47rp 


+  — (V  •  A)  =  -47rp 
c  dt 


(III.45) 


Equations  (III. 44)  and  (III. 45)  fully  describes  the  optical  field.  By  a  judicious  choice 
of  Coulomb  gauge,  V  •  A  =  0,  Equation  (III. 44)  becomes  the  full  wave  equation  in 
terms  of  transverse  current^° 

1  9^A  47r , 


V"A 


(III.46) 


dt‘^  c 

We  will  use  this  equation  to  determine  the  electron  beam  to  optical  mode  interaction. 

2.  FEL  Wave  Equation 

Eor  our  helical  undulator  given  by  Equation  (III. 3),  a  solution  to 

1  d‘^A  An  , 


V"A 


dt^ 


(III.47) 


^Vector  Identity  8,  [14] 
i°[10],  Chapter  6 
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can  be  written  as 


A(x,t) 


(111.48) 


where  E('x,t)  is  the  complex  optical  amplitude,  a  =  kz  —  ut  is  the  carrier  wave, 
i  =  (— z,  1,0)  is  optical  field’s  polarization  vector  for  circularly  polarized  fight,  and 
X  =  xi  +  yj  +  zk  is  the  position  vector.  Substituting  Equation  (III. 48)  into  (III. 47), 
we  get  the  following 


d‘^A  _  ee*“  d‘^E 
dx"^  k  dx"^ 
d‘^A  ee*“  d‘^E 


dy‘^  k  dy‘^ 


d^A 

ee*“ 

'd‘^E 

o-,  9E 

dz‘^ 

k 

_  dz"^ 

+  2ik— — 
dz 

-  k^E 

ee*“ 

'd'^E 

dE 

dt^ 

k 

dt^ 

1 

CC) 

3 

* 

1 

-  io^E 

(111.49) 

(111.50) 

(111.51) 

(111.52) 


For  a  laser  beam,  we  can  make  a  few  simplifications.  If  we  assume  that  optical  wave 
amplitude  and  phase  are  slowly  varying  along  the  2;  axis,  second  derivatives  in  2;  and  t 
can  be  dropped  compared  to  first  derivatives.  Substituting  these  back  into  Equation 
(III. 47)  and  using  u  =  kc,  we  arrive  at 


ee 


k 


_2  0-7  /  ^  1  ^ 

Vj_  +  2tk  (7^ — I —  — 


£  =  --3, 

c 


where  =  d‘^/dx‘^  +  d‘^/dy‘^  is  the  transverse  Laplacian  operator.  We  can  simplify 
the  expression  further  by  a  change  of  coordinates,  u  =  z  —  ct,  which  follows  the  light^^, 
and  recalling  dimensionless  time,  r  =  ct/L,  then 


A 

dz 

d 


dt 

,  d  Id 
so  that,  TT — \-  -tt 
oz  c  at 


du  d  dr  d  d 

dz  du  dz  dr  du 

dr  d  du  d  c  d  d 

dt  dr  dt  du  L  dr  ^du 

d  led  d 

du'^  c  L  dr  du 

L  dr 


(III.54) 


11 


Coordinate  transformation  is  also  known  as  method  of  characteristics.  See 


[15] 
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We  have  effectively  collapsed  2;  and  t  into  one  variable  r.  Multiplying  Equation 
(III. 53)  by  A:e*e*",  making  use  of  Equation  (111.54),  and  e*  •  e  =  2,  we  now  have 

= -27r-e-*“Jj.  •  r  (III.55) 

\LdT  J  \  c 

We  can  express  the  electron  current  density  as  a  sum  over  all  electrons 

Jx  =  — ec  E  /3_l5^^^(x  —  Xi(t))  (III. 56) 

i 

where  e  is  the  magnitude  of  electron’s  charge,  x;  is  the  electron’s  position  at  time  t, 
and  is  the  three-dimensional  delta-function.  Equation  (III. 11),  /3j.,  can  also  be 
expressed  as 

/3x  =  5R{-— (III.57) 

7 

where  5R  signifies  taking  the  real  part  of  the  complex  argument.  Combining  these 
equations.  Equation  (III. 55)  becomes 

/  1  f)  \ 

V‘\^2ik[-—]  E{r,  t)  =  -imkeKp  <  - >  (III.58) 

\  L  OT  J  \  7 

where  r  =  xi  +  yj.  We  have  made  use  of  Equation  (III. 14)  and 

p(x,t)=  f  —  Xi(t))dV  (III. 59) 

Jv  i 

The  average  <  ...  >  represents  an  average  over  sample  electrons  in  a  volume  element. 
If  we  multiply  Equation  (III. 58)  by  —AniNeKL"^ /'j'^mc^k,  recall  the  dimensionless 
optical  field  amplitude,  |a|  =  47rWeW|E'|L/7^mc^,  and  define  a  dimensionless  current 
density  j  =  87r^We^W^L^p/7^mc^,  we  have 

^  “(^’W)  =<  >  (III.60) 

when  7  7o,  a  reasonable  assumption  for  all  the  electrons  during  the  length  of  the 

interaction.  Looking  at  the  factors  multiplying  the  transverse  Laplacian  suggests  we 
define  dimensionless  transverse  coordinates  (where  the  tildes  refer  to  the  dimensionless 
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variables): 


(111.61) 
(111.62) 

We  can  now  recast  the  Laplacian  operator  into  a  dimensionless  Laplacian  operator, 
V^,  since 

92  _  k  I 

dx‘^  dy‘^  2L  dx"^  dy‘^ 

The  PEL  wave  equation  now  becomes, 

+  a(r,T)=<-je-*^  >  (III.64) 

where  a(f,T)  =  lale**^,  the  complex  dimensionless  optical  field  amplitude,  measures 
the  optical  field  strength  in  the  interaction  region;  the  dimensionless  current  j  mea¬ 
sures  the  interaction  between  the  electron  beam  and  the  optical  mode;  and  <  e  > 
measures  the  amount  of  electron  bunching.  If  j  is  large  (j  >>  tt),  the  optical  field  a 
changes  rapidly  so  we  have  high  gain.  If  j  is  small  (j  «  tt),  we  have  low  gain. 


3.  Gain 

Thus  far  in  our  development  of  EEL  physics,  we  have  studied  free  electrons 
interacting  with  laser  fight  through  a  Lorentz  force  interaction  in  the  undulator.  We 
have  described,  through  the  pendulum  equation,  the  process  of  physically  bunching 
the  electrons  in  phase  space  leading  to  coherent  radiation.  Einally,  we  have  coupled 
this  radiation,  through  Maxwell’s  equations,  with  the  electron  beam  co-propagating 
in  the  undulator  with  the  optical  beam  (see  Equation  (III.64)). 

We  must  now  discuss  the  modification  of  the  field  by  EEL  amplification.  In 
weak  fields,  |a|  <<  tt  and  low  gain  j  <<  tt,  we  can  apply  a  perturbation  expansion  of 
the  pendulum  equation  resulting  in  an  expression  for  optical  mode  amplication  [13] 
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where  the  filling  factor,  F,  is  the  ratio  of  areas  of  the  electron  beam  to  the  optical 
mode.  Gain  is  proportional  to  the  dimensionless  current  density,  j.  At  r  =  1,  the 
final  gain  is 


2  —  2  cos(z^o)  “  ^0  sin(z^o) 


(III.66) 


Equation  (III. 66)  is  shown  in  Figure  12.  Notice  the  anti-symmetric  nature  in 


Gain  vs.  v  ati  =  1 

0 


Figure  12.  Gain  curve  plotted  at  r  =  1. 


and  zero  gain  at  Vo  =  0.  The  gain  peaks  at  Pq  ~  2.6  with  peak  absorption  at 
Vo  ~  —2.6.  When  the  electrons  are  at  phase  velocities  slightly  above  resonance, 
the  amount  of  energy  transferred  from  the  electron  beam  to  the  optical  mode  is 
maximized  meaning  that  the  optical  mode  is  amplified  by  the  electron  beam.  In 
an  operating  FEL,  electrons  that  initially  go  through  the  undulator  spontaneously 
emit  in  a  spectrum  (Figure  13a)  centered  around  the  resonant  frequency.  The  weak 
field  low  gain  curve  (Figure  13b)  narrows  the  spectrum  over  many  passes  (Figure 
13c).  Over  many  passes,  the  laser  spectrum  narrows  due  to  mode  competition.  The 
top  of  Figure  13d  shows  a  narrow  linewidth  after  2500  passes  requiring  only  a  few 
microseconds  of  FEL  operation. 
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Figure  13.  FEL  coherence  as  it  evolves  from  a  broadband  spectrum,  (a)  Broad¬ 
band  spectrum  from  spontaneous  emission,  (b)  Low  field  gain  curve,  (c)  Narrowing 
broadband  spectrum  over  many  passes,  (d)  Final  power  spectrum  showing  narrow 
linewidth. 
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IV. 


OPTICAL  THEORY 


In  the  previous  section,  we  described  how  oscillating  electrons  co-propagating 
with  the  optical  field  gave  rise  to  amplified  laser  fight.  We  will  now  restrict  our 
discussion  to  the  the  laser  fight  itself  as  it  propagates  from  the  interaction  region  in  free 
space.  We  are  interested  in  the  transverse  optical  mode  structure  because  the  burn 
pattern  (oc  \a{x,  y)  p)  will  determine  the  intensity  delivered  to  a  beam  optic  or  a  target. 
Free  electron  lasers  (FEL)  primarily  operate  in  the  fundamental  Gaussian  mode,  but 
a  few  higher-order  modes  are  also  possible.  It  has  been  shown  in  simulations  and 
experiments  that  a  few  higher  order  modes  can  develop  when  the  FEL  is  operated 
with  mirror  shifts  and  tilts,  electron  beam  shifts  or  mirror  asymmetry.  Building  on 
previous  work  [5],  we  will  develop  the  mathematics  for  Hermite-Gaussian  modes  and 
a  means  of  finding  the  coefficients  for  a  given  waveform.  The  modal  analysis  described 
here  will  be  applied  to  EEL  mirror  distortions  which  is  being  investigated  for  the  first 
time. 


A.  OPTICAL  WAVE  EQUATION 


We  first  begin  our  discussion  of  optical  theory  by  describing  fight  propagation 
when  no  sources  are  present.  Recall  from  the  previous  chapter,  the  EEL  wave  equation 
in  the  undulator  is 


i  ~2  d 
+  — 


a(r,  r)  =<  —je  > 


4  ’  '  dT_ 

when  no  sources  are  present,  i.e.,  j  =  0  and  Equation  (IV. 1),  reduces  to 


(IV.l) 


.-9  ■ 

or 


a(f,  r)  =  0 


(IV.2) 


where  a  =  |a|e*'^,  the  complex  dimensionless  field  amplitude,  measures  the  optical 
field  strength.  In  terms  of  the  EEL  system,  the  following  description  is  applicable 
outside  of  the  undulator  or  when  the  electron  current  is  zero  inside  the  undulator. 
In  our  definitions  of  dimensionless  variables  x,  y,  and  r,  L  represented  the  length  of 
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the  undulator.  It  will  now  represent  the  length  that  the  optical  mode  travels.  Our 
task  now  will  be  to  find  solutions  of  Equation  (IV. 2)  to  investigate  transverse  mode 
properties. 

B.  FUNDAMENTAL  MODE 

A  form  for  the  solution  of  Equation  (IV. 2)  in  lowest  order  can  be  written  down 
as 

a(f,T)  =  aoe  (IV. 3) 

where  a(f ,  r)  is  the  dimensionless  complex  field  amplitude  proportional  to  the  electric 
field,  P  =  and  Gq  is  the  real  dimensionless  optical  amplitude  at  5;  =  0,y  = 

0,r  =  0.  Enactions  p(r)  and  qir)  are  yet  to  be  determined.  Rayleigh  length,  Zg  = 
Zg/L,  is  the  longitudinal  distance  from  the  mode  waist  at  which  the  optical  area 
doubles,  i.e.,  where  the  mode  radius  exceeds  the  value  at  the  waist  by  a  factor  of  a/2. 
It  is  a  measure  of  how  the  optical  mode  expands  (see  Eigure  14).  Equation  (IV. 3) 


a 


b 


0x1 

Eigure  14.  Cross  section  of  a  Gaussian  mode  propagating  from  r  =  0  — >■  1.  The  mode 
with  a  smaller  Rayleigh  length  expands  more  quickly,  (a)  Zg  =  0.20,  =  0.0  (b) 

Zg  =  0.8,  Tw  =  0.0 
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is  the  dimensionless  analogue  of  a  form  used  by  [16],  et  al.  Since  we  are  solving  a 
second-order  differential  equation,  we  must  specify  initial  conditions  to  determine  a 
unique  solution.  We  stipulate  that  at  r  =  0,  a(r,  0)  =  age  ,  an  initially  transverse 
Gaussian  wavefront.  The  initial  conditions  then  require: 

p(0)  =  0 

g(0)  =  1  (IV.4) 

Our  job  will  be  to  find  p(t)  and  qir). 

We  now  take  derivatives  of  Equation  (IV. 3)  resulting  in: 


d'^c 


-2 


=  ^  + 


dx‘^  \Zoq  zlq^ 


(IV.6) 


d‘^a  _  /  “2  ^  Ay 


dy‘^ 


Zoq  zlq^ 


(IV.6) 


da  (  °  x"^  q  y"^  q\ 

(Recall  (.?.)  =  d{. . .) / dr .)  Combining  Equations  (IV. 5),  (IV.6),  (IV. 7)  into  (IV. 2) 
and  dividing  through  by  a,  we  are  arrive  at 


2  ,  2  ^  Ay^  ^  °  ^ 

- h  ~  - h  „  -\-  Al  \  —  P  ^  +  T — 5-  1  —  0 

Zoq  z^q^  Zoq  z^q^  \  Zoq^  Zoq^ 


Making  use  of  =  x"^  +  y"^,  we  have 

.2 


Zoq^  \Zo 


—  +  iq  ]  =  \  —  +  ip 


Zoq 


(1V.8) 


(IV.9) 


Eor  this  equation  to  hold  for  any  at  some  particular  time  r,  both  terms  in  paren¬ 
theses  must  be  identically  zero.  Solving  the  left  hand  side  for  q(T) 


Q  =  — 


i 

Zo 


qir)  = 


-dr' 

Zo 


IT 

=  qoA-  — 
Zo 


(IV.IO) 
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Applying  the  initial  conditions,  Equation  (IV. 4),  we  find 


q{T)  =  I  +  — 
^0 


(IV.ll) 


Solving  the  right  hand  side  of  Equation  (IV. 9)  for  pij)  and  substituting  Equation 
(IV.ll),  we  find  P{t)  to  be 

(IV.12) 


p= 


+ 


T 


Zl  +  zl  + 

Integrating  with  respect  to  r  and  applying  initial  conditions,  we  have 


(IV.13) 


j-  I  ^2  I  2;2 

p{t)  =  i  arctan  —  +  In  '  ° 


(IV.14) 


Substituting  Equations  (IV.ll)  and  (IV. 14)  back  into  Equation  (IV. 3),  we  find 


a(r,  t)  =  ttge 


-i  arctan  In 


(IV.15) 


We  can  clean  up  this  mess  if  we  make  a  few  substitutions.  Erom  Equation  (IV.ll), 
we  have 


I 

q{^) 


I  +  V 

Zo 


T 


—  l- 


Zq  'V  ~ 

°  Zo 


Zq  ~ 

°  Zo 


With  the  benefit  of  hindsight,  we  can  make  the  following  definition: 


w‘^{t)  =  Zo+^ 

^0 


(IV.16) 


(IV.17) 


Equation  (IV.  17)  represents  the  optical  mode  radius  as  shown  in  Eigure  15.  The  term 
Zo  is  the  dimensionless  analogue  of  the  dimesional  Rayleigh  length.  The  derivation 
has  assumed  that  the  beam  waist  occurs  at  r  =  0,  that  is,  t(;(0)  =  Wo  =  But 
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the  beam  waist  may  be  positioned  anywhere  along  r.  It  is  a  trivial  matter  to  make 
a  time  translation  in  the  equations  such  that,  r  — >  (r  —  r^).  With  Equation  (IV. 17), 
we  can  re-express  1/q  in  terms  of  w  and  Wg-  Equation  (IV. 16)  now  becomes 


1 

Qi^) 


T 


—  I 


(IV.18) 
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Eigure  15.  Cross  section  of  a  Gaussian  mode  as  it  propagates  from  r  =  0  — >■  1 
illustrating  the  optical  mode  radius,  =  0.3,  Zg  =  0.3 


Returning  to  Equation  (IV.  15)  and  substituting  Equation  (IV.18),  the  optical 


field  becomes^ 


Wg _ 

a(r,T)  =  ag-——re 
w(t) 


„-2  . 
%- 


-i  arctan  ^ 


(IV.19) 


Einally,  we  can  write  down  the  fundamental  mode  of  the  optical  field  as 


Amplitude 


a(r,r)  =  ttg 


Wr, 


-e  *  e 


-arctan  — 
•20 


W[T] 


Phase 


(IV.20) 


^The  square  root  term  can  be  expressed  more  compactly  as  Wojw 
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where  the  mode  waist  is  at  =  0.  For  any  arbitrary  location  of  the  mode  waist,  we 


have 


Amplitude 


a(r,  T)  =  tto 


Wn 


WiT  — 


r  (t-tw) 

u^w‘^('r-'rw) 


arctan  - 


(IV.21) 


Phase 


As  can  be  seen,  the  amplitude  falls  off  as  and  is  clearly  Gaussian  while  the 

phase  term  includes  both  f  and  r  dependence.  Equations  (IV. 20)  and  (IV.21)  repre¬ 
sents  the  lowest-order  transverse  electromagnetic  mode.  Figure  16  shows  the  surface 


0  0  0  0 


Figure  16.  Real  and  imaginary  parts  of  complex  Gaussian  mode  at  r  =  0  with 
Tyj  =  0.5. 

plots  of  the  real  and  imaginary  parts  of  Equation  (IV.21).  The  real  and  imaginary 
parts  of  a(f,  r)  provides  us  with  the  amplitude  and  phase  of  the  optical  field  since 

|a(f,r)|  =  y/5R{a}2  +  (IV.22) 

and 

t)  =  arctan  (IV.23) 

Ja|a| 

where  (f)  is  the  optical  phase.  Burn  patterns  are  recognizable  through  the  power  in  the 
optical  field,  i.e.,  the  optical  field  amplitude  squared.  In  the  following  discussions  we 
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will  only  plot  the  amplitude  of  the  mode  instead  of  the  power  (or  amplitude  squared) 
because  we  would  lose  detail  in  the  plots. 


initial  surface  plot  final  surface  plot 


0  0  y  0  0 


Figure  17.  Initial  (r  =  0)  and  final  (r  =  1)  surface  amplitudes  for  the  lowest  order 
mode  where  yellow  represents  highest  largest  amplitude. 


X  0  X  1  X 


Figure  18.  Three  panels  shows  cross  sections  of  a  Gaussian  mode  as  it  propagates 
from  r  =  0  — >■  1.  Zo  =  0.25  =  0. 


Figure  17  is  a  surface  plot  for  a  Gaussian  mode  at  r  =  0  and  r  =  1.  Figure 
18  is  composed  of  three  panels.  The  first  panel  shows  a  cross  section  {x  vs.  y)  of  the 
Gaussian  mode  at  r  =  0.  The  middle  panel  shows  a  cross  section  {y  vs.  r)  as  the 
mode  propagates.  The  last  panel  shows  the  cross  section  of  the  mode  at  r  =  1  (x  vs. 
y).  Notice  how  a  small  Rayleigh  length  has  spread  the  mode  considerably. 
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C.  HIGHER  ORDER  MODES 


In  order  to  account  for  more  complicated  variations  along  transverse  directions, 
we  can  generalize  the  trial  solution  Equation  (IV. 3): 


a(r,Tj  =  ttog 


X 


W[T 


h 


W[T 


-{p{t)  + 


loq{-r) 


(IV.24) 


The  functions  g  and  h  will  allow  for  a  more  intricate  variance  along  transverse  direc- 
tions.  Again,  we  demand  the  same  initial  conditions,  namely,  a(f ,  0)  =  age  .  This 
implies,  p(0)  =  0,  g(0)  =  1,  and  g{x/wo)  =  h{y/wo)  =  1  since  the  initial  Gaussian 
form  already  expressed  in  the  factor  As  before,  our  job  now  will  be  to  find  the 

unknown  functions,  p{t),  qir),  g{x/w{T)),  and  h{y/w{T))  at  later  times.  But  before 
we  begin  taking  derivatives,  let  us  consider  derivatives  of  g  (and  by  analogy  h)  to 
illustrate  the  chain  rule: 


and 


dg  dg  d  f  x\ 

dx  d{^)  dx  \tc/ 


1  dg 


d‘^g  _  d  dgi 
dx‘^  dx  dx 


dg  dg  d  /  X  \ 

dr  d{^)  dr  \w) 

_  ~  dg  /  1  \  Stc 


(IV.25) 


(IV.26) 


(IV.27) 


As  you  can  see,  the  primes  over  g  denote  a  derivative  with  respect  to  its  argument. 
This  little  fact  will  become  handy  later.  Taking  derivatives  of  Equation  (IV.24)  results 


in: 
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e 


(IV.28) 


d‘^a 

dx"^ 


ig'hx 


WZoQ 


2gh  Aghx‘^ 
ZoQ  zlq^ 


-(p(t)  +  - 


Zq{-r)  > 


d‘^a 

dy"^ 


igh'y 


WZoQ 


2gh  4:ghy‘^ 
ZoQ  zlq^ 


-{p{t)+ 


ioq{-r) 


(IV.29) 


1  °  ^ 

-^g'h  w  X 


1  °  _  ° 

-^gh'  w  y-ghp  -[ 


ghx"^ 

ZoQ^ 


Q  4 


ghy"^ 

ZoQ^ 


Q 


e 


(p(t)+ 


2 

^og(T) 


(IV.30) 


Substituting  Equations  (IV.28),  (IV.29),  and  (IV.30)  into  (IV. 2),  and  dividing 
through  by  ttoghe  we  have  the  following  unpleasant  monster: 


Ig  Ax  g  w  ~g  1  h 

- TT - Ai—x—  +  — ^ - 

g  WZoQ  g 

^  ^  's. 


Ay  h'  w  ^h' 

~2-  •  ~2  U  + 

w‘^  g  w‘^  n  WZoQ  n  n 


X  only 


y  only 


+  r" 


+ 


Ai 

ZoQ^ 


- - AiP  =  0  (IV.31) 

ZoQ 

" - V - ' 


only 


no  r  dependence 


In  order  for  Equation  (IV.31)  to  have  a  solution  at  a  given  time  r  (i.e.  a  transverse 
plane),  the  x  dependent  terms  should  equal  a  constant,  let’s  say,  —a.  The  y  dependent 
terms  should  equal  another  constant,  let’s  say  — /5.  The  dependent  terms  should 
be  identically  zero  since  it  must  be  valid  for  at  any  location  in  the  transverse  plane. 
Einally,  the  terms  that  don’t  depend  on  r,  should  equal  a  +  /5,  so  that  the  left  hand 
side  sums  to  zero.  Solving  the  r  equation,  we  have 


+ 


Ai 


ZoQ^ 

4f^  r  1  ° 

—  +  iQ 


ZoQ^  '-Zo 


=  0 

=  0 


(IV.32) 

(IV.33) 


Eor  this  equation  to  hold  for  any  r,  the  term  in  bracket  must  be  identically  zero.  But 
we  have  seen  this  requirement  before  in  the  fundamental  mode  derivation.  Applying 
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the  initial  conditions,  we  find  that 


q{T)  =  1  +  — 
^0 


(IV.34) 


Returning  to  the  x  dependence  of  Equation  (IV. 31),  we  multiply  through  by 
gw"^  and  collect  terms  in  derivatives  of  g,  giving  us 


g”  —  Ax  (  - \-  i  w\  g'  +  aw^g  =  0 

\Zoq  J 


(IV.35) 


Let  us  evaluate  the  parenthetical  term,  recalling  Equations  (IV. 17),  (IV. 18),  and 

O 

noting  that  w=  r/wlw,  we  find 


So  now  we  have 


w  . 

- \-  1  w  \ 

w 

(?i- 

Zoq  J 

wi 

1 

w 

// 

9 

Ax  , 

- —g  + aw  g 

(IV.36) 


(IV.37) 


The  above  equation  is  close  to  the  form  of  the  Hermite  equation  [12],  y"  —  2xy'  +  \y  = 
0.  Recall  that  the  primes  denote  a  derivative  with  respect  to  the  argument  of  the 
function.  So  a  judicious  change  in  variables  will  allow  us  to  recast  Equation  (IV.37) 
in  precisely  the  form  of  Hermite’s  equation.  Let 


\/2x  _ 

^  ^ - y  X  =  —j=. 

w  1/2 


=  V2d{^] 

w 


Using  Equations  (IV. 38)  and  (IV. 39)  on  (IV.37),  results  in 

9g^  ^dg  aw‘^ 

a|  -  2«af  +  — « =  “ 


(IV.38) 

(IV.39) 

(IV.40) 


A  power  series  solution  will  solve  Equation  (IV.40)  uniquely  with  initial  conditions. 
But  for  the  series  to  terminate,  we  must  require  aw"^ /2  to  be  equal  to  a  nonnegative 
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even  integer.  With  this  requirement,  the  power  series  solution  terminates  into  the 
well  known  Hermite  polynomials: 


9(0  = 


where  m 


0,1, 2, 3... 


(IV.41) 


There  will  be  a  similar  solution  for  the  y  terms  of  Equation  (IV. 31),  i.e.,  h{r])  = 
Hn{ri)  =  Hn{V^y/w)  where  n  =  0,1,  2,  3...  The  first  ten  Hermite  polynomials  are 
fisted  in  Table  II. 


Ho{x)  =  1 

Hi{x)  =  2x 

H^lx)  =  Ax"^- 2 

H‘^{x)  =  %x^  —  12x 
H4^(x)  =  16x‘^  —  +  12 

H^{x)  =  ?>2x^  —  IGOa:^  +  120a: 

H^{x)  =  64a:®  -  480a:4  +  720a:2  -  120 
H^\x)  =  128a:^  -  1344a:®  +  3360a:3  -  1680a: 

H^{x)  =  256a:®  -  3584a:®  +  13440a:4  -  13440a:2  +  1680 

H^{x)  =  GWa:*^  -  9216a:^  +  48384a:®  -  80640a:®  +  30240a: 


Table  II.  Hermite  Polynomials  from  [12] 


Returning  to  Equation  (IV. 31),  we  are  left  with  finding  a  solution  for  pir) 
Eor  Hermite  solutions,  a  and  /5  are  constrained 


aw^  „  4m 

=  2m  — >■  q;  =  — r 


2 


^  .  4n 

=  2n  — >■  /5  = 


2 

o 

Solving  for  P,  with  the  constraints,  and  Equations  (IV. 34)  and  (IV. 17), 


-4:ip=  a  +  13 

ZoQ 


(IV.42) 


P  =  +  n)  +  ^ 

ZoQ 

=  i{m  +  n  +  1)  ° 


+ 


T 


K  + 


(IV.43) 
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Integrating  and  applying  initial  conditions  for  p 


p{t) 


7" 

i{m  +  n  +  1)  arctan  —  +  In 
^0 


(IV.44) 


So  we  now  have  all  the  pieces  to  Equation  (IV. 24).  Putting  it  altogether,  Equations 
(IV. 34),  (IV. 41),  and  (IV.44)  into  (IV. 24),  and  after  a  little  tidying  up,  we  arrive 
finally  at  an  expression  for  the  (m,n)  mode  for  the  laser  optical  field: 


Amplitude 


a^,Ji,T)  =  ao'^H^  (  ^  (  ^  ]  e-& 

W  \w  \w  - - - 

/  \  y  Phase 


(IV.45) 


where  the  phase,  (f)m,n,  is  given  by 


(r,T)  = 


T 


—  (m  +  n  +  1)  arctan  — 


(IV.46) 


Eigure  19  shows  the  “burn  patterns”  for  various  modes. 


Eigure  19.  |aoo|,  |aoi|,  |cf3o|,  l^ail  (Left  to  right;  Top  to  bottom) 
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D.  NORMALIZATION 


If  we  attempt  to  superpose  modes  with  high  values  of  m  or  n,  the  Hermite 
polynomials  with  Gaussian  envelopes  for  these  higher  modes  have  large  coefficients 
as  shown  in  Table  II.  To  place  all  the  modes  on  an  equal  footing  we  will  normalize 
each  of  them  to  the  power  in  the  fundamental  which  is  given  by 


-Po,o  — 


00  noo 


ttQfi'^dxdy 


00  ^  —00 


Substituting  Equation  (IV. 20)  we  have 


P,.o  - 


2  /'oo 


9  POO  .2 

/  2y^ 


e 


dx  e 

J  —OO 


(IV.47) 


(IV.48) 


so  then 

n.o  =  (IV.49) 

The  power  in  an  arbitrary  mode  is 


P  = 

m,n 


00  POO 


^m,n  dxdy 


(IV.50) 


oo  — oo 


Substituting  Equation  (IV. 45)  we  have 


P  = 

4  m,n 


f'OO  ^2  ^,2 

al-z^H^Hle^^dxdy 

-oo  ^ 


(IV.51) 


To  evaluate  this  integral,  we  make  use  of  the  fact  that  Hermite  polynomials  are 
orthogonal  on  the  range(— oo,oo)  with  respect  to  a  weighting  function  [12], 


(IV.52) 


where  the  Sij  is  the  Kronecker  delta.  To  properly  make  use  of  the  orthogonality 
condition.  Equation  (IV.52),  we  must  first  massage  Equation  (IV.51)  into  the  correct 
form.  Consider  the  factor 


(IV.53) 
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(IV.54) 


If  we  make  the  following  change  of  variables,  let 

\/2x  _  aw 


a  —  — ; —  — )■  X  —  — ;= 
w 


so  that, 


dx  =  —;=da 


(IV.55) 


Then  Equation  (IV. 53)  becomes 

~  poo 

^  /  Hm{afe-^"da  (IV.56) 

V  2  J  —  oo 

Employing  the  orthogonality  condition.  Equation  (IV. 52)  when  m'  =  m,  we  have 

J  oo 

and  similarly  for  the  y  factors.  The  power  in  the  (m,n)  mode  then  becomes 

Pm,n  =  |a>22™2"m!n!  (IV.58) 

=  Po,o2™2”m!n!  (IV.59) 


H„ 


fQ;)^e 


-^"da  =  w.  -2^m\ 


(IV.57) 


So  the  power  in  the  (m,n)  mode  can  be  normalized  to  the  power  in  the  fundamental 
by  dividing  it  through  by  2"*2”m!n!.  A  normalized  expression  for  the  optical  field  in 
the  (m,n)  mode  is  then 


4  —  _ , 

r- — - \ — r’ 

V2™2”m!n! 


(IV.60) 


This  expression  places  all  the  higher-order  modes  on  an  equal  footing;  all  modes  have 
the  same  power  as  the  fundamental. 


E.  FINDING  COEFFICIENTS 

The  total  optical  field  is  a  superposition  of  all  the  separate  normalized  com¬ 
ponents,  i.e., 

OO 

Ar(f,T)  =  ^C'^,n^m,n(f,T)  (IV.61) 

m,n 

where  C,fn^n  is  the  normalized  coefficient  for  the  (m,n)  mode.  Each  term,  or  com¬ 
ponent,  is  a  solution  to  the  wave  equation.  So,  the  sum  of  terms  is  also  a  solution. 
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Moreover,  the  Hermite-Gaussian  modes  form  a  complete  orthogonal  basis  set  [16]  so 
that  any  laser  field  satisfying  the  wave  equation  can  be  accurately  represented  by 
Equation  (IV. 61)  with  appropriate  coefficients  Cm,n-  Our  task  now  is  to  determine 
the  unknown  normalized  coefficients  Cm,n-  As  with  a  Fourier  series  decomposition 

of  a  function,  we  can  employ  the  Hermite  analog  of  Fourier’s  Trick^.  To  determine 

* 

the  coefficients  for  an  optical  field.  At,  we  multiply  Equation  (IV. 61)  by  Am',n'  and 
integrate  over  all  x  and  y  giving  us 

/OO  poo  ^  poo  poo  ^ 

I  At  Ajjii dxdy  /  /  ^  ]  Cjji^nAm,n  A^t dxdy  (IV. 62) 

OO  j —OO  J—ooJ—oo^^ 

where  the  asterisk  denotes  taking  the  complex  conjugate.  Expanding  out  the  right 
hand  side,  we  have 


OO  ~ 

Er  ^  ^ 

V2™2”m!n!  w 


^  ^Orr  rr  - 

V2™2”m!n!  w 


(IV.63) 

Notice  the  term  that  involves  the  product  of  the  phases  from  the  unprimed  and  primed 
optical  field.  Recalling  Equation  (IV. 46),  we  then  have 


_  —  m)  +  (n'—  n))arctanj^] 


(IV.64) 


This  equation  does  not  depend  on  x  or  y,  so  we  can  pull  it  out  of  the  integral.  Re¬ 
turning  to  Equation  (IV.63)  and  employing  the  orthogonality  of  Hermite  polynomials 
when  m'  =  m  and  n'  =  n,  we  have 


r>00  poo 


At  dxdy 


-OO  ^  — OO 


'2"*2”m!n! 


rAyj'^ 


fi;W|2"n! 


(IV.65) 


Combining  this  result  with  the  left-hand  side  of  Equation  (IV. 62),  we  have 


r*  — 

'^m,n 


OO  poo 


TTClgWg  J —OO  ■-> —OO 


At  A^„  dxdy 


(IV.66) 


^[14],  Chapter  3 
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This  equation  is  an  expression  for  the  normalized  coefficients  for  our  optical  field  At- 
Our  next  task  will  be  to  write  computer  code  to  determine  the  coefficients  of  any 
waveform. 
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V. 


SIMULATIONS 


In  this  chapter  we  will  describe  how  numerical  simulation  can  solve  the  non¬ 
linear  differential  equations  that  have  been  developed.  Non-linear  differential  equa¬ 
tions  are  exceedingly  difficult  to  solve  analytically  but  are  relatively  “easy”  with  a 
computer. 


A.  MODELING  THE  FREE  ELECTRON  LASER 

In  order  to  model  the  PEL  we  must  simultaneously  solve  two  coupled  differ¬ 
ential  equations:  the  pendulum  equation,  Equation  (III. 32), 

OO 

c=  |a|  cos(C  +  (/>) 


and  the  EEL  wave  equation.  Equation  (III. 64), 


i  ~ 


d 


a(r,T)  =<  -je 


-iC 


> 


The  pendulum  equation  updates  electron  phases  and  phase  velocities  and  the  wave 
equation  propagates  the  light  self-consistently  with  the  electron  beam  interaction. 
The  Naval  Postgraduate  School’s  Directed  Energy  and  Electric  Weapons  Center  uses 
an  Apple  Xserve  cluster  to  model  the  EEL  (see  Eigure  20).  The  cluster  is  specifically 
designed  to  model  the  EEL  in  four  dimensions  as  well  as  model  EEL  performance  due 
to  perturbations  or  test  validity  of  new  designs. 

In  order  to  numerically  solve  Equation  (111.64),  fast  Eourier  transforms  are 
used  in  the  code  [17].  An  output  for  a  recent  simulation  is  given  in  Eigure  21. 
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Figure  20.  NFS  Apple  Xserve  cluster  with  64  nodes  /  128  processors. 
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nx=800,  nt=BO<  iip=2S000,  tniir=19,  1^=20  Fri  May  26  IB; 33; 46  2006 

FEL  WAVEFRONT  EVOLUTICN 


|a(x,y, 10)1 


j=60.0 

Vo=12.B 

<j  =0.22 

a^=0.22 
a  =0.31 
a,J=0.31 

1  ®p=1.2 

cr  =0.26 

r  e”  =0 

yo 

a;=10 

VO.B 

z  =0.42 

rj^=0.0 

r“=40 

1  e=3e-07 

r  =9.B 

«n=9 

N=30 

-0.B00%(c) 

P(n) 


l.ft<10^ 


-n/2  4  3x72 


[ST  (h)  0-31 

liuimifM 


147 

|a(x,y,-9)|  2 

(0,y,-9) 

y 

t 

147 

(i) 

(i) 

(0,y,10) 


-147  X  147-3  3-147  x  147  -12 

P=l. 666403  (2%),  G=0.14B  (97%),  r|=0.0144  (4%),  av=28.B 
Ginn=0.0279,  Giinx=0.311,  Gf=0.14B 
F=0.31,  G(th)=2,  G(ext)=O.OB8,  A7/y=0.07B7 
v0=0.81  (0.6B),  vl=14  (IS),  v2=14  (IB) 

M^=0.9B,  c^(2,0)=0.31,  c^(0,0)=0.31,  c^(4,0)=0.27 
f4iins=67,  ebins=32 


12 


|c(»n,ii)| 


Figure  21.  FEL  simulation  output,  (a)  Optical  mode  evolution  over  500  passes  at 
right  mirror  (b)  Surface  plot  of  the  optical  field  at  the  right  mirror  after  last  pass; 
(c)  User  defined  parameters  (see  Table  (HI));  (d)  Cross  section  of  optical  mode  and 
electron  beam  as  it  propagates  from  r  =  0  — >■  1  on  last  pass;  (e)  Evolution  of  electron 
beam  distribution  in  phase  velocity  over  500  passes;  (f)  Electron  phase  space  at  r  =  1 
after  last  pass;  (g)  Optical  power  evolution  over  500  passes;  (h)  Optical  gain  evolution 
over  500  passes;  (i)  Burn  pattern  of  optical  mode  at  left  mirror  after  last  pass;  (j) 
Optical  phase  for  at  left  mirror  after  last  pass;  (k)  Burn  pattern  of  optical  mode  at 
right  mirror  after  last  pass;  (1)  Optical  phase  for  at  right  mirror  for  last  pass;  (m) 
Hermite-Gaussian  coefficients  for  right  mirror  (Code  is  developed  in  this  chapter) 


49 


Symbol 

Name 

Equation 

j 

Electron  current  density 

j  =  87r^We^iL^L^p/7^mc^ 

Electron  phase  velocity 

^0  =  [{k  +  ko)(5z  -  k]L 

ujp 

Betatron  frequency 

ujp  =  2ttNK/'j 

Electron  beam  waist  position 

Tp  =  Zp/L 

Vo 

Electron  beam  shift 

Vo  =  YosJt^/LX 

dyo 

Electron  beam  tilt 

Oyo  =  Q\/L'k/X 

a 

Initial  optical  field  amplitude 

\a\  =  47rNeK\E\L/'j‘^mc^ 

Zo 

Rayleigh  length 

zq  Zq!  l 

Dimensionless  mirror  radius 

e 

Mirror  edge  loss 

See  [18] 

Qn 

Cavity  quality  factor 

Qn  =  l/(loss  per  pass) 

Electron  beam  radius 

(7x  Xg-\^Tl  J  LX 

^0X1  ^0y 

Electron  beam  angular  spread 

O’l'e 

Spread  in  electron  phase  velocity 

due  to  emittance 

See  [18] 

Spread  in  electron  phase  velocity 

due  to  energy  spread 

See  [18] 

At! 

Optical  mode  waist  position 

At!  Yyj  !  L 

Th 

Right  mirror  hole  radius 

Th  =  Rh'jT^I^L 

Tc 

Mirror  radius  of  curvature 

> 

II 

Table  III.  User-defined  PEL  simulation  parameters 


Each  node  of  the  Apple  cluster  can  run  a  simulation  such  as  the  one  shown 
above  at  a  different  values  of  Pq  simultaneously.  We  can  then  determine  gain  and 
extraction  values  at  the  optimum  value  of  Pq-  This  simulates  what  the  EEL  would  do 
by  mode  competition. 

The  simulation  also  provides  us  with  an  array  describing  the  optical  field, 
a{x,y).  As  mentioned  in  the  previous  chapter,  the  optical  field  is  a  superposition  of 
Hermite-Gaussian  basis  set  with  appropriately  chosen  coefficients.  Mirror  distortion 
or  electron  beam  shifts  and  tilts  can  cause  the  optical  field  to  develop  higher-order 
modes.  We  will  now  analyze  these  modes  by  finding  their  coefficients  numerically. 
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B.  FINDING  HERMITE-GAUSSIAN  COEFFICIENTS  NU¬ 
MERICALLY 

Equation  (IV. 66)  is  the  analytical  expression  for  the  normalized  Hermite- 

Gaussian  (HG)  coefficients.  We  need  to  recast  this  equation  in  a  form  that  is 

* 

computer-friendly.  Since  At  and  Am,n  are  complex,  we  can  write: 


At  =  Ap  +  iAp 

(V.i) 

A*  —4’’  —  1  4* 

^m.n  ^^m.n 

(V.2) 

where  Alj.  and  Alp  are  the  real  and  imaginary  parts  of  the  total  normalized  optical 
field  and  A^^^^  and  A^^^^  are  the  real  and  imaginary  parts  of  the  normalized  (m,n) 
mode.  Equation  (IV. 66)  becomes 

2 


r  — 

.n 


(V.3) 


Expanding  out  the  factors  and  separating  the  coefficients  into  their  real  and  imaginary 
parts,  we  have 


and 


Cl 


Cl 


2 

7ra>2 

2 

TralWg 


+  -^T  Kt.n)  rfi 

-  At  A'„)  dx  dy 


(V.4) 

(V.6) 


+  (CiJ-  (V.6) 

So,  let  us  first  describe  the  procedure  for  developing  a  numerical  modal  analysis  tool 
in  the  G  programming  language: 


1.  Gonstruct  Hermite-Gaussian  basis  arrays. 

2.  Read  in  or  form  a  user-defined  optical  field. 

3.  Determine  HG  coefficients  via  Equations  (V.4),  (V.5),  and  (V.6) 
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Our  expression  for  the  optical  field,  Equation  (IV. 61)  is  a  sum  of  the  normalized  basis 
at  some  time  r.  Recall  that  the  optical  field  is  a  complex  number  so  we  must  keep 
track  of  the  real  and  imaginary  parts.  These  analytical  expressions  are  continuous  in 
X  and  y  so  we  recast  the  expression  for  the  (m,n)  mode  into  their  real  and  imaginary 
parts  over  the  discrete  sites  in  the  array.  Therefore,  in  pseudo-code  we  have 


A''[m][n][xj\[yk]  = 


n 

W 


n 

W 


V2x. 


H„ 


V2yk\  Jh 

e  *2  cos 


w 


w 


—  {m  +  n  +  1)  arctan  ^  )  (V.7) 


A^[m][n][xj][yk]  = 


y/2Xj 


w 


w 


Jh 

e 


sm 


—  {m-\-  n  -\-l)  arctan  —  )  (V.8) 


Zo 


where  and  A^  are  four-dimensional  arrays  indexed  by  the  mode  number  and  posi¬ 
tion,  given  by  (m,n)  and  (7j,  yk).,  respectively.  We  calculate  the  basis  set  for  the  first 
several  modes  and  store  these  values  in  their  arrays.  We  then  read  in  the  optical  field 
that  is  passed  as  an  array  from  the  main  EEL  simulation  or  form  a  user-defined  field^ 
by  a  linear  combination  of  normalized  Hermite-Gaussian  modes  using  the  arrays  that 
been  calculated. 

We  determine  the  HG  coefficients  via  Equations  (V.4),  (V.5),  and  (V.6).  We 
approximate  the  double  integral  by  nested  loops  over  m,  n,  and  over  the  number  of 
array  sites,  nx.  In  pseudo-code, 

^  nx  nx 

C"[m][n]  =  Y  {^T[xj\[yk]  A^[m][n][xj\[yk]  + 

"  "  k  j 

(V.9) 

and 


nx  nx 


G*[m][n]  =  —  [xj][yk]  A''[m][n][xj][yk]  -  AJj,[xj][yk]  A^[m][n][xj][yk]) 

fv.ioi 


7ra2fi;2  ^ 


user-defined  field  can  be  used  as  a  check  on  the  code  to  ensure  that  it  is  operating  properly. 
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The  square  of  the  magnitude  of  the  HG  coefficient  for  the  (m,n)  mode  is  then  given 
by  Equation  (V.6),  that  is, 

C[m][n\^  =  iC^[m\[n]f  +  iC^[m][n\f  (V.ll) 


It  is  convenient  to  normalize  the  square  of  the  magnitudes  to  the  power  of  the 
total  field  At-  The  power  in  the  total  field  is  given  by 


Pt  = 


poo 

8 

f  —00  ^ 

—  00 

(V.12) 


Substituting  Equation  (IV. 61)  and  working  out  the  integrals  we  find  that 

00 

=  (V.13) 

m^n 

Recall  that  the  power  in  the  fundamental  mode  is  given  by  Po,o  =  so  that 

the  total  power  in  the  optical  field  is  then  (in  terms  of  the  fundamental) 


Pt  =  -Po.0  5] 

m^n 

The  fractional  power  spectrum,  is  then 


S[m][n  = 


G[m][n]^ 

WM 


(V.14) 


(V.16) 


where  jC[i][jY  is  the  sum  of  all  the  squares  of  the  normalized  mode  coefficients. 

To  illustrate  how  the  code  works,  consider  the  following  examples  illustrated  in 
Eigures  22  and  23.  In  Eigure  22  (a),  the  diagnostic  correctly  identifies  the  fractional 
power  spectrum  as  50%  in  both  Aq^2  and  ^2,0  modes.  Incidentally,  the  superposition 
of  equally  weighted  Aq^2  and  ^2,0?  with  no  relative  phase  difference  between  them, 
results  in  a  pure  L?  Laguerre-Gaussian  mode^. 

In  Eigure  22  (b),  the  diagnostic  correctly  identifies  the  fractional  power  spec¬ 
trum  for  the  superposition  of  Aq^,  ^1,1,  and  A3  3  with  Go,o  =  1,  Gi^i  =  2,  and  ^3^3  =  3 


^To  see  this:  ask  Professor  Colson  and  he  will  buy  you  a  beer. 
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as  7%,  28%,  and  64%,  respectively.  A  screen  capture  of  Cmn.c  is  shown  in  Figure 
23  for  this  case.  As  you  can  see,  the  program  outputs  the  Rayleigh  length  and  mode 
waist.  It  then  displays  the  mode  radius  at  r  =  0,  the  power  in  the  fundamental,  and 
total  power  in  the  arbitrary  field.  The  fractional  power  spectrum  is  then  calculated 
and  displayed.  The  code  then  propagates  the  optical  field  in  free-space  and  performs 
the  same  analysis  on  the  final  waveform. 

In  Figure  22  (c),  the  diagnostic  correctly  identifies  the  fractional  power  spec¬ 
trum  for  the  superposition  of  equally  weighted  Aq^,  Ao,2,  and  Ao,4  modes.  We  will  see 
a  similar  combinations  of  even  ‘m’  and  ‘n’  valued  modes  when  we  investigate  mirror 
distortions. 

The  diagnostic  developed  for  this  thesis  for  Hermite-Gaussian  modal  analysis 
was  then  incorporated  into  the  main  NFS  FEL  simulation  code  as  a  function  call. 
In  the  next  chapter,  we  will  see  this  diagnostic  in  action  when  we  simulate  mirror 
distortions. 
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Code  results 

50%  Aq2 
50%  A20 


7%Aoo 
29%  All 
64%  A33 


33%  Aqo 
33%  Ao2 
33%  Aq4 


Figure  22.  Graphics  from  Cmn.c  output  was  drawn  with  MatLab  (a)  Superposition 
of  equally  weighted  4o,2  and  42,0-  (b)  Superposition  of  4o,o,  ^1,1,  and  ^3^3  with 
Go,o  =  1,  =  2,  and  6*3^3  =  3.  (c)  Superposition  of  equally  weighted  Aq^,  ^0,2, 

and  4o,4. 
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Files  opened 

Initializations 

zo  =  0.5  tanw  =0.5 

Loading  np  basis  sets  at  tan  =  0. 

w(0)  =  1 

P_00(tan  =  0)  =  78.5398 
P_1(0)  =  1099.56 

Ontpntting  initial  wavefront  to  file 
Calculating  initial  Cinns 

NoriTBlized  initial  cirm  amplitudes  squared 
Cmn2[0][0]  =  0.0714286 
Cmn2[l][l]  =  0.285714 
Cmn2[3][3]  =  0.642857 

Propagating  waveform. . . 

w(l)  =  1 

P_00(1)  =  78.5398 
P_I(1)  =  1119.46 

NbriTBlized  initial  cmn  amplitudes  squared 
Cmn2[0][0]  =  0.0702047 
Cmn2[l][l]  =  0.28182 
Cmn2[3][3]  =  0.639894 


Closing  files... 

sjianl  05-250:  ~/Desktop/Cnrrent  Code  rvigilpbg4$ 

I 


Figure  23.  Output  from  diagnostic  tool  Cmn.c 
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VI. 


MIRROR  DISTORTIONS 


A.  COLD  CAVITY  THEORY 

In  the  pursuit  of  a  weapons  class  FEL,  high  power  and  high  intensity  can  lead 
to  mirror  distortions  if  heat  deposited  by  the  laser  light  to  the  mirror  substrate  isn’t 
removed  by  cooling  in  a  uniform  manner.  Distortions  then  cause  higher-order  modes 
to  appear.  It  is  important  to  quantify  the  fractional  power  spectrum  among  the 
higher-order  modes  because  doing  so  uniquely  identifies  the  laser  field  and  may  allow 
for  corrective  optics  to  single  out  the  best  modes  for  FEL  lethality.  We  will  investigate 
three  different  cases  of  mirror  distortions  on  the  outcoupling  mirror:  hyperbolic, 
ellipsoidal,  and  spherical.  Hyperbolic  and  ellipsoidal  distortions  can  be  classified  as 
a  mirror  astigmatism.  It  occurs  when  the  radius  of  curvature  of  the  spherical  mirror 
along  an  axis  differs  by  some  small  amount  when  compared  to  the  perpendicular  axis. 
For  simplicity  we  chose  these  axes  to  fine  up  with  x  and  y.  Spherical  distortion  is 
not  astigmatism  because  both  axis  are  changed  in  exactly  the  same  amount  and  in 
the  same  direction.  Spherical  distortions  do  change  the  Rayleigh  length.  Before  we 
consider  the  effects  of  mirror  distortions  on  FEL  operation,  let  us  explore  cold-cavity 
effects  (i.e.,  no  electron  beam  present)  of  mirror  “errors”  to  get  a  flavor  for  how  the 
FEL  might  respond.  Keep  in  mind  that  cold-cavity  theory  discussed  in  this  section 
keeps  both  mirrors  the  same.  The  dimensionless  radius  of  curvature  of  the  two  ideal 
mirrors  determines  the  Rayleigh  length,  Zo  [18],  and  is  given  by^ 


"  2  ^  Trr, 


(VI.l) 


where  is  the  dimensionless  mirror  separation.  Factoring  out  r^/S  and  defining  a 
small  parameter  k  =  45^ /r^,  we  then  have 


=  y(i  +  «) 


(VI.2) 


^Symmetric  confocal  cavity.  See  [19],  chapter  12. 
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where  k  <<  1  for  >>  Zg.  Typically,  ~  10  and  Zg  ~  0.1,  so  k  4  x  10^"^.  What 
if  an  actual  mirror  differs  from  the  ideal  mirror  by  a  small  amount?  This  could  occur 
because  mirror  heating  has  distorted  both  mirrors,  for  instance.  We  can  quantify  this 
distortion  as  follows 

fg  fc(l  +  e)  (VI. 3) 

where  e  is  a  small  deviation  from  the  original  radius  of  curvature  fg-  Combining  these 
relations  we  have 

f,  =  ^(l  +  K)(l  +  e)  (VI.4) 

fg  =  ^(1  +  K  +  e  +  Ke)  (VI. 5) 

Since  k  and  e  are  small,  we  can  approximate  this  equation  to  first  order  giving  us 

fg  ~  -^(1  +  K  +  e)  (VI. 6) 

We  now  ask,  “What  is  going  on  with  the  mode  radius  at  the  waist  and  at  the 
outcoupling  mirror  with  respect  to  k  and  r^?”  The  waist  is  given  by 


We  can  solve  Equation  (VI.  1)  for  Zg  giving  us 


and  substitute  that  into  our  expression  for  the  waist  giving  us 


Substituting  Equation  (VI. 6),  we  now  have 


(  Z'm  \  2  /  \  i 

Wg  ^  (K  +  e)4 


(VL7) 

(VL8) 

(VL9) 

(VLIO) 

(VLll) 
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Let  us  now  consider  the  mode  radius  at  the  mirror.  The  mode  radius  is  given 
by  Equation  (IV. 17).  At  the  outcoupling  mirror,  r  =  r^/S,  we  then  have 


w{Tm/2)  =  W, 


(VI.12) 


Substituting  Equations  (VI. 6)  and  (VI. 8)  in  for  Zo  and  after  a  little  algebra  we  arrive 
at 


Wm  =  \[^  (VI.  13) 

V  K  +  e 

In  both  Equations  (VI. 11)  and  (VI. 13),  there  is  a  problem  when  the  mirror  distortion 
e  — >■  —K.  Let  us  plot  these  two  equations  with  respect  to  e/n.  Erom  Eigure  24  we 


Waist  radius 


e/K 


Mode  radius  at  mirror 


(a) 


(b) 


Eigure  24.  Cold-cavity,  (a)  Waist  radius  variation  (b)  Mode  radius  variation  at  the 
outcoupling  mirror. 


can  see  that  if  the  error,  e  is  less  than  —k  the  waist  radius  and  mode  radius  becomes 
imaginary.  If  e  =  —  k,  the  waist  radius  becomes  zero  and  the  mode  radius  at  the 
mirror  becomes  infinity.  So,  as  far  as  cold-cavity  analysis  predicts,  mode  radius  can 
be  sensitive  to  mirror  distortions  that  are  on  the  order  of  —k.  So  we  may  expect 
actual  EEL  operation  to  follow  a  similar  trend. 
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We  chose  to  simulate  the  effects  of  mirror  distortions  on  two  different  JLab 
oscillator  configurations:  Electromagnetic  (EM)  wiggler  and  STI  Optronics  wiggler 
designs.  The  following  table  highlights  their  major  differences. 


EM  wiggler  STI  wiggler 


Undulator  period,  Xg  (cm) 

8 

5.5 

Undulator  parameter,  K 

0.5 

1.36 

Undulator  length,  L  (m) 

2.3 

1.65 

No.  Periods,  N 

29 

30 

Mirror  Separation,  S  (m) 

30 

50 

Quality  factor, 

25 

9 

Rayleigh  length,  Zg  (cm) 

90 

70 

K  (%) 

0.35 

0.2 

Table  IV.  JLab  EM  and  STI  wiggler  comparison 


B.  HYPERBOLIC  DISTORTION 

In  the  case  of  hyperbolic  distortion,  we  vary  the  mirror  distortion  on  the 
outcoupling  mirror  in  the  following  way: 

Afcy  =  -Afcx  (VI.  14) 

where  Afcy  and  Afcx  are  the  changes  of  outcoupling  mirror  radius-of-curvature  in  the 
X  and  y  directions,  respectively.  Eirst  let  us  consider  the  EM  wiggler  design.  The 
results  for  gain  and  extraction  are  shown  in  Eigure  25.  Notice  how  both  the  gain  and 
extraction  graphs  are  symmetric  about  zero  e  as  should  be  expected  for  a  hyperbolic 
distortion.  Although,  cold  cavity  analysis  predicts  catastraphic  results  when  mirror 
distortions  are  less  than  —k,  we  see  that  with  EEL  operation  both  the  gain  and 
extraction  curves  continues  well  beyond  this  cold  cavity  limit  (— k  is  indicated  by  the 
red  tick  mark). 

Notice  also  a  peak  in  gain  and  extraction  as  the  magnitude  of  e  is  increased. 
It  seems  counterintuitive  that  increasing  hyperbolic  mirror  distortion  would  improve 
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Figure  25.  (a)  Gain  trend  for  hyperbolic  distortion  for  the  JLab  EM  wiggler.  (b) 
Extraction  trend  for  hyperbolic  distortion  for  the  JLab  EM  wiggler 


EEL  operation.  This  can  be  explained  if  we  look  at  at  Eigure  25  (b).  If  we  compare 
the  burn  patterns  for  e  =  0  and  e  =  0.5%,  we  can  see  that  the  mode  shape  at  the 
mirror  (yellow)  is  wider  for  the  e  =  0.5%  case  than  for  the  e  =  0%  case.  If  the 
mode  is  more  spread  out  at  the  mirror,  we  can  conclude  that  the  mode  was  more 
tightly  focused  in  the  undulator  leading  to  a  greater  interaction  between  the  optical 
mode  and  the  electron  beam  leading  to  greater  gain  and  extraction.  This  beneficial 
effect  does  not  go  on  indefinitely.  Eventually  greater  distortion  will  cause  the  mode 
to  become  so  large  that  the  mode  will  be  clipped  by  the  outcoupling  mirror  reducing 
gain  and  extraction.  An  example  of  a  mode  being  clipped  is  shown  in  Eigure  25  (a) 
for  the  e  =  —  2%  case.  We  also  show  the  modal  coefficients  that  comprise  this  data 
point.  Notice  that  the  ‘n’  modes  are  weighted  more  heavily  than  the  ‘m’  modes  since 
the  burn  pattern  is  spread  out  in  the  y  direction  and  that  only  even  modes  are  present 
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since  the  original  wave  shape  is  assumed  to  be  a  fundamental  Gaussian.  Only  the 
even  modes  overlap  the  origin  (a:  =  0,  y  =  0)  where  the  electron  beam  is  amplifying 
the  light. 


Figure  26.  (a)  Gain  trend  for  hyperbolic  distortion  for  the  JLab  STI  wiggler.  (b) 
Extraction  trend  for  hyperbolic  distortion  for  the  JLab  STI  wiggler 

We  next  considered  the  STI  wiggler  design.  The  results  for  gain  and  extraction 
are  shown  in  Figure  26.  As  with  the  EM  wiggler  case,  we  see  that  the  gain  and  the 
extraction  are  symmetric.  The  gain  is  much  larger  in  the  STI  wiggler  design  due  to  a 
larger  value  of  the  undulator  parameter,  K.  The  mirror  separation  is  greater,  thereby 
yielding  a  somewhat  smaller  value  for  k  {kem  =  0.35%  vs.  ksti  =  0.2%).  We  also  see 
a  rise  in  gain  and  extraction  for  distortions  on  the  order  of  |k|  as  in  the  EM  wiggler 
case.  Notice  in  Figure  26  (b)  for  e  =  0.5%,  how  the  coefficients  are  strongly  weighted 
in  the  ‘m’  modes  due  to  a  spreading  out  in  the  x  direction.  Again,  only  the  even 
modes  are  driven  by  the  electron  beam  at  the  origin. 
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C.  ELLIPSOIDAL  DISTORTION 


In  ellipsoidal  distortion,  we  vary  the  mirror  error  on  the  outcoupling  mirror  in 
the  following  way: 

vary  Af^  ,  Af^x  =  0  (VI.  15) 

As  before,  let  us  consider  the  EM  wiggler  results.  Notice  in  Figure  27,  both  the  gain 


Figure  27.  (a)  Gain  trend  for  ellipsoidal  distortion  for  the  JLab  EM  wiggler.  (b) 
Extraction  trend  for  ellipsoidal  distortion  for  the  JLab  EM  wiggler 


and  extraction  curves  are  not  symmetric,  as  should  be  expected  for  an  asymmetric 
distortion.  Both  gain  and  extraction  curves  have  the  same  general  shape  as  Figure 
24  (a).  Unlike  the  cold  cavity  case,  the  EEL  operates  well  beyond  e  <  —k,.  As 
e  is  decreased  further,  higher-order  modes  appear  as  shown  in  the  coefficient  chart. 
Conversely,  as  e  is  increased  (e  >  —k),  the  mode  remains  near  the  fundamental  mode. 

In  Figure  28,  we  have  the  results  for  ellipsoidal  distortions  in  the  STI  wiggler. 
Both  curves  are  once  again  asymmetric.  Notice  how  the  gain  curve  reaches  its  max- 
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'0  ^  ^  S 

t  (%)  (=  A  r^frj  (b) 


Figure  28.  (a)  Gain  trend  for  ellipsoidal  distortion  for  the  JLab  STI  wiggler.  (b) 
Extraction  trend  for  ellipsoidal  distortion  for  the  JLab  STI  wiggler 


imum  at  e  =  —1.5%  which  is  several  times  greater  than  cold-cavity  theory  predicts 
the  resonator  should  fail.  The  extraction  curve  is  not  as  smooth  as  the  EM  wiggler 
extraction  curve.  But  the  general  shape  is  still  present,  i.e.,  maximum  nigh  — e  and  a 
relatively  flat  extraction  in  the  positive  e  direction. 

D.  SPHERICAL  DISTORTION 

In  spherical  distortion,  we  vary  the  mirror  distortion  on  the  outcoupling  mirror 
in  the  following  way: 

vary  Afcy  =  Af^x  (VI.  16) 

Spherical  distortion  is  not  a  misnomer  even  though  the  mirror  remains  cylindrically 
symmetric,  because  the  highly  reflecting  mirror  is  not  altered.  The  high  power  laser 
beam  is  assumed  to  have  caused  a  symmetric  distortion,  as  in  Equation  (VI. 16).  Res- 
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onator  cavity  instabilities  are  introduced  as  the  distortion  is  increased  or  decreased. 
Notice  in  Figure  29,  the  burn  patterns  are  cylindrically  symmetric.  Again,  we  see 
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Figure  29.  (a)  Gain  trend  for  spherical  distortion  for  the  JLab  STI  wiggler.  (b) 
Extraction  trend  for  spherical  distortion  for  the  JLab  STI  wiggler 


good  gain  and  extraction  well  into  the  regime  where  cold-cavity  theory  predicts  a 
non-functioning  FEL.  Modal  decomposition  reveals  an  equal  weighting  of  lower-order 
modes  yielding  burn  patterns  with  central  holes  in  them.  Spherical  distortion  varies 
the  effective  cavity  Rayleigh  length.  In  Eigure  30,  we  notice  the  burn  patterns  caused 
by  varying  Rayleigh  length  are  cylindrically  symmetric  and  comprised  of  equally 
weighted  modes  similar  to  spherical  distortions.  In  Eigure  30  (a),  gain  increases  as 
we  go  to  shorter  Rayleigh  lengths  because  the  mode  is  more  focused  in  the  interac¬ 
tion  region  (see  Eigure  14)  increasing  the  value  of  the  filling  factor.  Since,  gain  is 
proportional  to  the  filling  factor,  gain  goes  up.  Shorter  Rayleigh  lengths  expands  the 
mode  on  the  mirrors.  In  Eigure  30  (b),  the  green  curve  shows  a  decreasing  intensity 
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Figure  30.  (a)  Gain  trend  for  Rayleigh  length  variation  in  the  JLab  STI  wiggler.  (b) 
Extraction  and  mirror  intensity  trends  for  Rayleigh  length  variation  in  the  JLab  STI 
wiggler 


with  decreasing  Rayleigh  length.  This  is  good  for  EEL  optics  because  the  decreased 
intensity  will  mitigate  mirror  distortions  caused  by  uneven  heating  or  cooling. 
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VII. 


CONCLUSION 


In  this  thesis,  we  have  reviewed  free  electron  laser  theory  and  for  the  first 
time  developed  analytical  solutions  to  quantify  Hermite-Gaussian  higher-order  modes, 
developed  a  diagnostic  for  modal  analysis,  and  investigated  tolerance  limits  on  mirror 
distortions.  Hermite-Gaussian  modes  form  a  complete  and  orthogonal  basis  set  that 
can  fully  describe  any  optical  field  and  is  a  natural  basis  set  when  describing  mirror 
distortions  that  do  not  possess  cylindrical  symmetry. 

Modal  analysis  is  a  useful  tool  for  describing  mode  quality.  We  are  able  to 
identify  the  fractional  power  spectrum  in  higher-order  modes.  This  information  may 
be  used  by  an  adaptive  optics  system  to  single  out  the  best  modes  for  propagation 
and  lethality. 

Tolerance  limits  on  mirror  distortions  are  several  times  greater  than  what  cold 
cavity  theory  predicts.  This  is  due  to  the  electron  beam  not  allowing  the  optical  mode 
to  shrink  to  zero  size.  The  electron  beam  stabilizes  the  operation  of  the  PEL.  Sim¬ 
ulations  indicate  JLab  EEL  should  tolerate  mirror  distortions  several  percent,  which 
is  much  greater  than  cold-cavity  predictions.  Einally,  we  recommend  the  following 
future  work: 

•  Investigate  Laguerre-Gaussian  modes  as  another  complete  and  orthogonal  ba¬ 
sis  set  to  include  azimuthal  dependence. 

•  Incorporate  both  Hermite-Gaussian  and  Laguerre-Gaussian  modal  decompo¬ 
sition  in  the  complete  four  dimensional  EEL  simulation  program. 

•  Utilizing  the  diagnostic  tool,  investigate  higher-order  modes  that  appear  due 
to  mirror  shifts/tilts  and  electron  beam  shifts/tilts. 

•  Ghoose  a  catchy  and  snappy  name  for  NPS  EEL  simulation  software.... something 
other  than  wavens. 


67 


THIS  PAGE  INTENTIONALLY  LETT  BLANK 


68 


APPENDIX  A.  CODES 


We  present  the  following  programs  used  for  finding  Hermite-Gaussian  coef¬ 
ficients  of  an  arbitrary  optical  field  and  the  MatLab  scripts  used  in  plotting.  Line 
numbers  are  included  for  convenience. 


1.  CMN.C 

/  *  start  of  Cmn.c  */ 

1  //  LCDR  Ricardo  Vigil 

2  //  Thesis  work 

3  //  June  2006  graduation  date 

4  //  13  major  revisions 

5 

6  /*  Program  functions: 

7  -  Combines  up  to  ten  user  defined  Hermite  Gaussian  modes 

8  to  form  an  initial  optical  field 

9  -  Decomposes  the  initial  optical  field  into  its  constituent  HG  modes 

10  -  Propagates  the  waveform  in  free  space 

11  -  Decomposes  the  propagated  waveform  into  its  constituent  HG  modes. 

12 

13  input  files:  Cmn.in 

14 

15  output  files: 

16  initial. out  -  (x,y)  initial  amplitude  in  optical  field,  i.e,  at  tau  =  0 

17  final. out  -  (x,y)  final  amplitude  in  optical  field,  i.e.,  at  tau  =  1 

18  ytau.out  -  (y,t)  amplitude  cross  section  in  (y,  tau  =  0  to  1) 

19 

20  Matlab  graphics  file:  Coeffs.m 

21  Creates  5  graphics 

22  1 .  Amplitude  burn  pattern  at  tau  =  0 

23  2.  Amplitude  surface  plot  at  tau  =  0 

24  3.  Amplitude  burn  pattern  at  tau  =  1 

25  4.  Amplitude  surface  plot  at  tau  =  1 

26  5.  Propagation  cross  section  in  y  vs  tau 

27 

28  */ 

29 

30  #include  <stdio.h> 

31  #include  <stdlib.h> 
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32  #include  <math.h> 

33 

34 

35  /*  these  arrays  are  defined  outside  mainO  due  to  segmentation  faults  */ 

36 

37  double  Abinitialr[10] [10] [200] [200] ;  //  real:  (m,n)  tau  =  0 

38  double  Abinitiali [10] [10] [200] [200] ;  //  img:  (m,n)  tau  =  0 

39  double  Abfinalr[10] [10] [200] [200] ;  //  real:  (m,n)  tau  =  1 

40  double  Abf inali [10] [10] [200] [200] ;  //  img:  (m,n)  tau  =  1 

41 

42  int  main(){ 

43 

44  double  P0=0.0,  P1=0.0,  rsq=0.0,  dx=0.0,  dy=0.0,  C=0.0,  wi=0.0,  wf=0.0; 

45  double  Ao=0.0,  zo=0.0,  wo=0.0,  tau=0.0,  tauw=0.0,  dtau=0.0,  phi=0.0; 

46  double  c2sumlnitial=0 . 0,  c2sumFinal=0 . 0,  templ=0.0,  temp2=0.0; 

47  double  Hx=0.0,  Hy=0.0,  ampinitial=0 . 0 ,  ampf inal=0 . 0 ,  ampy=0.0,  K=0.0; 

48  double  pi=3 . 1415926 ,  u=0.0,  v=0.0; 

49 

50  double  cinitialr[10] [10] ;  //  real:  inital  coeff  for  the  mth  nth  mode 

51  double  cinitiali [10] [10] ;  //  img:  intial  coeff  for  the  mth  nth  mode 

52  double  c2maginitial [10] [10] ; //  initial  magnitude  squared  for  (m,n)  mode 

53  double  cfinalr[10] [10] ;  //  real:  final  coeff  for  the  mth  nth  mode 

54  double  cf inali [10] [10] ;  //  img:  final  coeff  for  the  mth  nth  mode 

55  double  c2magf inal [10] [10] ;  //  final  magnitude  squared  for  (m,n)  mode 

56  double  iS[10][10];  //  initial  power  spectrum  ratio  for  (m,n)  mode 

57  double  fS[10][10];  //  final  power  spectrum  ratio  for  (m,n)  mode 

58 

59  double  Powerinitial [10]  [10] ;  //Power  in  (m,n)  normalized. 

60  double  Powerf inal [10] [10] ;  //Power  in  (m,n)  normalized. 

61 

62  double  Atr [200]  [200] ;  //real:  total  real  optical  field 

63  double  Ati  [200]  [200] ;  //img:  total  img  optical  field 

64  double  Atr_old[200] [200] ;  //old  values:  real 

65  double  Ati_old[200] [200] ;  //old  values:  img 

66 

67  int  i=0, j=0,k=0,l=0,nx=0,¥=0,x=0,y=0; 

68  double  c0=0.0,cl=0.0,c2=0.0,c3=0.0,c4=0.0; 

69  double  c5=0.0,c6=0.0,c7=0.0,c8=0.0,c9=0.0; 

70  int  m0=0,ml=0,m2=0,m3=0,m4=0,m5=0,m6=0,m7=0,m8=0,m9=0; 

7 1  int  n0=0 , nl=0 , n2=0 , n3=0 , n4=0 , n5=0 , n6=0 , n7=0 , n8=0 , n9=0 , count 1=0 ; 

72  unsigned  int  f; 

73  unsigned  int  factorial (unsigned  int  a); 
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75  int  m[10]  ,  n[10]  ; 

76  double  c[10]; 

77 

78  FILE  *input,  ^initial,  *final,  *ytau; 

79  input  =  fopenC'Cmn.  in"  ,  "r")  ; 

80  initial  =  fopenC'initial.out"  ,  "w") ; 

81  final  =  f openC'f  inal .  out"  ,  "w") ; 

82  ytau  =  f  openC'ytau.  out"  ,  "w") ; 

83 

84  printf ("Files  opened\n") ; 

85 

86  fscanf  (input ,  "“/olf  Zd  "/d"  ,&c0,&m0,&n0) ; 

87  fscanf  (input ,  "“/olf  “/od  “/d"  ,&cl,&ml,&nl)  ; 

88  fscanf  (input ,  "“/olf  “/d  “/d"  ,&c2,&m2,&n2)  ; 

89  fscanf  (input ,  "“/olf  “/d  “/d"  ,&c3,&m3,&n3) ; 

90  fscanf  (input ,  "“/olf  “/d  “/d"  ,&c4,&m4,&n4)  ; 

91  fscanf  (input ,  "“/olf  “/d  “/d"  ,&c5,&m5,&n5) ; 

92  fscanf  (input ,  "“/olf  “/d  “/d"  ,&c6,&m6,&n6)  ; 

93  fscanf  (input ,  "“/olf  “/d  “/d"  ,&c7,&m7,&n7) ; 

94  fscanf  (input ,  "“/olf  “/d  “/d"  ,&c8,&m8,&n8)  ; 

95  fscanf  (input ,  "“/olf  “/d  “/d"  ,&c9,&m9,&n9) ; 

96  fscanf  (input ,  "“/olf  “/d  “/d"  ,&Ao,&W,&nx) ; 

97  fscanf  (input ,  "“/olf  “/If  “/If  "  ,&zo  ,&tauw,&dtau) ; 

98 

99  /*  initializations - */ 

100 

101  printf ("Initializations\n") ; 

102 

103  c[0]=c0;  c[l]=cl;  c[2]=c2;  c[3]=c3;  c[4]=c4; 

104  c[5]=c5;  c[6]=c6;  c[7]=c7;  c[8]=c8;  c[9]=c9; 

105 

106  m[0]=m0;  m[l]=ml;  m[2]=m2;  m[3]=m3;  m[4]=m4; 

107  m[5]=m5;  m[6]=m6;  m[7]=m7;  m[8]=m8;  m[9]=m9; 

108 

109  n[0]=n0;  n[l]=nl;  n[2]=n2;  n[3]=n3;  n[4]=n4; 

110  n[5]=n5;  n[6]=n6;  n[7]=n7;  n[8]=n8;  n[9]=n9; 

111 

112  dx  =  (1.0*¥)/(1.0*nx) ; 

113  C  =  dtau/(4.0*dx*dx) ; 

114  wo  =  sqrt(zo); 

115  printf  ("zo  =  “/gXt  tauw  =  “/gXn" ,  zo,tauw); 
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116 

117  for(k=0;k<10;k++) 

118  for(l=0;l<10;l++)  { 

119  Powerinitial [k]  [1] =0 . 0 ; 

120  Powerfinal[k] [1]=0.0; 

121  cinitialr[k] [1]  =  0.0; 

122  cinitiali [k] [1]  =  0.0; 

123  c2maginitial[k] [1]  =  0.0; 

124  cfinalr[k] [1]  =  0.0; 

125  cf inali [k] [1]  =  0.0; 

126  c2magfinal[k] [1]  =  0.0; 

127  iS[k]  [1]  =  0.0; 

128  fS[k]  [1]  =0.0;  } 

129 

130  f or (x=0 ; x<nx; X++) 

131  for(y=0;y<nx;y++)  { 

132  Atr[x]  [y]  =  0.0; 

133  Ati  [x]  [y]  =  0.0; 

134  Atr_old[x] [y]  =  0.0; 

135  Ati_old[x] [y]  =  0.0;  } 

136 

137  for(k=0;k<10;k++) 

138  for(l=0;l<10;l++) 

139  for (i=0; i<nx; i++) 

140  for(j=0; j<nx; j++){ 

141  Abinitialr[k]  [1]  [i]  [j]  =  0.0; 

142  Abinitiali [k] [1] [i] [j]  =  0.0; 

143  Abfinalr[k]  [1]  [i]  [j]  =  0.0; 

144  Abf inali [k]  [1]  [i]  [j]  =  0.0;  } 

145 


146  /*  Initializaitons  complete  - */ 

147  /*  Loading  basis  set - */ 

148 

149  printf ("\nLoading  up  basis  sets  at  tau  =  0.\n\n") 


150 

151  tau  =  0.0; 

152  wi  =  sqrt(  zo  +  (tau-tauw) * (tau-tauw) /zo  ); 

153  printf  ("w(y„g)  =  y„g\n"  ,tau,  wi)  ; 

154 

155  for(k=0;k<10;k++) 

156  for(l=0;l<10;l++) 

157  for(x=0;x<nx;x++) 
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158  for(y=0;y<nx;y++){ 

159  u  =  (sqrt(2)/wi)*(x  -  nx/2  +  0.5)  *  dx; 

160  V  =  (sqrt(2)/wi)*(y  -  nx/2  +  0.5)  *  dx; 

161  if  (k==0)  Hx=1.0; 

162  if  (k==l)  Hx=2*u; 

163  if  (k==2)  Hx=4*pow(u,2)-2; 

164  if  (k==3)  Hx=8*pow(u,3)-12*u; 

165  if  (k==4)  Hx=16*pow(u,4)-48*pow(u,2)+12; 

166  if  (k==5)  Hx=32*pow(u,5)-160*pow(u,3)+120*u; 

167  if  (k==6)  Hx=64*pow(u,6)-480*pow(u,4)+720*pow(u,2)-120; 

168  if  (k==7)  Hx=128*pow(u,7)-1344*pow(u,5)+3360*pow(u,3)-1680*u; 

169  if  (k==8)  Hx=256*pow(u,8)-3584*pow(u,6)+13440*pow(u,4)-13440* 

170  pow(u,2)+1680; 

171  if  (k==9)  Hx=512*pow(u,9)-9216*pow(u,7)+48384*pow(u,5)-80640* 

172  pow(u,3)+30240*u; 

173  if  (1==0)  Hy=1.0; 

174  if  (1==1)  Hy=2*v; 

175  if  (1==2)  Hy=4*pow(v,2)-2; 

176  if  (1==3)  Hy=8*pow(v,3)-12*v; 

177  if  (1==4)  Hy=16*pow(v,4)-48*pow(v,2)+12; 

178  if  (1==5)  Hy=32*pow(v,5)-160*pow(v,3)+120*v; 

179  if  (1==6)  Hy=64*pow(v,6)-480*pow(v,4)+720*pow(v,2)-120; 

180  if  (1==7)  Hy=128*pow(v,7)-1344*pow(v,5)+3360*pow(v,3)-1680*v; 

181  if  (1==8)  Hy=256*pow(v,8)-3584*pow(v,6)+13440*pow(v,4)-13440* 

182  pow(v,2)+1680; 

183  if  (1==9)  Hy=512*pow(v,9)-9216*pow(v,7)+48384*pow(v,5)-80640* 

184  pow(v,3)+30240*v; 

185 

186  rsq  =  ((x-nx/2+0.5)*(x-nx/2+0.5)+(y-nx/2+0.5)*(y-nx/2+0.5))*dx*dx; 

187  phi  =  (rsq* (tau-tauw) / (wo*wo*wi*wi) - (k+1+1) *atan( (tau-tauw) /zo) ) ; 

188 

189  Abinitialr [k] [1] [x] [y]=Ao*wo*Hx*Hy*exp(-rsq/(wi*wi))*cos(phi)/wi* 

190  sqrt (1/ (pow (2 , k) *pow (2 , 1) *f actorial (k) *f actorial (1) ) ) ; 

191 

192  Abinitiali [k] [1] [x] [y] =Ao*wo*Hx*Hy*exp(-rsq/ (wi*wi) ) *sin(phi)/wi* 

193  sqrt (1/ (pow (2 , k) *pow (2 , 1) *f actorial (k) *f actorial (1) ) ) ; 

194 

195  Powerinitial[k] [1]  =  Powerinitial [k] [1] + 

196  Abinitialr [k] [1] [x] [y] *Abinitialr [k]  [1] [x] [y]  + 

197  Abinitiali [k] [1] [x] [y] *Abinitiali  [k]  [1]  [x]  [y] ;  } 

198 

199  /*  This  serves  as  a  check  to  ensure  modes  add  up  in  the  expected  way  */ 
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200 

201  printf  ("P_00(tau  =  0)  =  °/og\n"  ,Powerinitial  [0]  [0] ) ; 

202 

203 

204 

205  /*  forming  initial  At  @  tan  =  0 - */ 

206  for(i=0; i<10; i++) 

207  for(x=0;x<nx;x++) 

208  for(y=0;y<nx;y++)  { 

209  Atr[x][y]  =  Atr[x]  [y]  +  c  [i]  *Abinitialr [m[i] ]  [n[i]  ]  [x]  [y]  ; 

210  Ati[x][y]  =  Ati[x][y]  +  c  [i]  *Abinitiali  [m[i]  ]  [n[i]  ]  [x]  [y]  ;  } 

211 

212  /*  finding  total  power  in  At  */ 

213 

214  for(x=0;x<nx;x++) 

215  for(y=0;y<nx;y++) 

216  PO  =  PO  +  Atr[x]  [y]  *Atr[x]  [y]  +  Ati  [x]  [y]  *Ati  [x]  [y]  ; 

217 

218  printf  ("P_T(0)  =  y.g\n",P0); 

219  printf ("Outputting  initial  wavefront  to  file\n\n"); 

220 

221  /*  Output  At  cross  section  to  file  and  saving  old  values*/ 

222 

223  f or (x=0 ; x<nx; X++)  { 

224  for(y=0;y<nx;y++)  { 

225  ampinitial  =  sqrt(  Atr  [x]  [y]  *Atr  [x]  [y]  +  Ati  [x]  [y]  *Ati  [x]  [y]  ); 

226  Atr_old[x]  [y]  =  Atr[x]  [y]  ; 

227  Ati_old[x]  [y]  =  Ati[x][y]; 

228  fprintf  (initial,  "“/If  ",  ampinitial)  ;  } 

229  fprintf (initial,  "\n") ;  } 

230 

231  /*  Cmn  calculation  for  initial  waveform - */ 

232 

233  printf ("Calculating  initial  Cmns\n"); 

234 

235  K  =  (2*dx*dx)/(pi*zo) ; 

236  for(k=0;k<10;k++) 

237  for(l=0;l<10;l++) 

238  for(x=0;x<nx;x++) 

239  for(y=0;y<nx;y++)  { 

240  cinitialrEk]  [1]  =  cinitialr[k]  [1]  +  K  * 

241  (Atr[x]  [y]*Abinitialr[k]  [1]  [x]  [y]  + 
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242  Ati [x] [y] *Abinitiali [k] [1] [x]  [y] ) ; 

243 

244  cinitiali [k] [1]  =  cinitiali [k]  [1]  +  K  * 

245  (Ati[x]  [y]*Abinitialr[k]  [1]  [x]  [y]  - 

246  Atr [x] [y] *Abinitiali [k] [1] [x]  [y] ) ; } 

247 

248 

249  /*  calculating  Cmn  squared  */ 

250  for(k=0;k<10;k++) 

251  for(l=0;l<10;l++)  { 

252  c2maginitial[k] [1]  =  (cinitialr[k] [1] *cinitialr  [k]  [1]  + 

253  cinitiali [k] [1] ^cinitiali [k] [1] ) ; 

254  c2suinlnitial  =  c2sumlnitial  +  c2maginitial  [k]  [1]  ;  } 

255 

256  printf ("\nNormalized  initial  cmn  amplitudes  squared\n"); 

257 

258  for(k=0;k<10;k++) 

259  for(l=0;l<10;l++)  { 

260  iS[k][l]  =  c2maginitial[k]  [l]/c2sumlnitial; 

261  if  (iS[k][l]  >  0.01)  printf  ("Cmn2 [y„d]  ["/od]  =  y„g\n"  ,k,l,iS[k]  [1] )  ;  } 

262 

263 

264  /*  propagation  of  At  from  tau  =  0  to  tau=  1  */ 

265 

266  printf ("\nPropagating  waveform. . .\n\n") ; 

267 

268  for(tau=0.0;tau<=0.9999;tau=tau+dtau)  {/*  start  tau  loop  */ 

269 

270  /*  wave  equation  over  x,y  and  tau  */ 

271  for(y=l;y<nx-l;y++) 

272  for(x=l;x<nx-l;x++){ 

273  Atr[x]  [y]  =  Atr  [x]  [y] -C*  (Ati_old[x+l]  [y] +Ati_old[x-l]  [y]  + 

274  Ati_old[x] [y+1] +Ati_old[x] [y-l]-4.0*Ati_old[x] [y] ) ; 

275 

276  Ati[x][y]  =  Ati[x]  [y]+C*(Atr_old[x+l]  [y]+Atr_old[x-l]  [y]  + 

277  Atr_old[x] [y+l]+Atr_old[x] [y-l]-4.0*Atr_old[x]  [y]) ;  } 

278 

279  /*  saving  old  values  */ 

280  for(x=0;x<nx;x++) 

281  for(y=0;y<nx;y++){ 

282  Atr_old[x]  [y]  =  Atr[x]  [y]  ; 

283  Ati_old[x]  [y]  =  Ati[x][y];  } 
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284 

285  /*  y  vs .  tau  output  */ 

286  if  (countl==100)  { 

287  count 1=0; 

288  for(y=0;y<nx;y++)  { 

289  ampy  =  sqrt  (Atr [nx/2] [y] *Atr [nx/2]  [y]  +  Ati [nx/2] [y] *Ati [nx/2] [y] ) ; 

290  fprintf  (ytau,  "“/og  ",ampy);  }  } 

291 

292  fprintf (ytau," \n") ; 

293 

294  countl  =  countl  +1;  }  /*end  tau  loop  */ 

295 

296  /*  output  for  final  waveform  */ 

297  for(x=0;x<nx;x++)  { 

298  for(y=0;y<nx;y++){ 

299  ampfinal  =  sqrt(  Atr[x]  [y]  *Atr  [x]  [y]  +  Ati  [x]  [y]  *Ati  [x]  [y]  ); 

300  fprintf  (final ,  "“/og  ",  ampfinal);  } 

301 

302  fprintf (final, "\n") ;  } 

303 

304 

305  /*  Calculating  cmns  for  final  waveform  (i.e.,  after  propagation) -  */ 

306  /*  Must  create  basis  set  at  tau  =  1.0  */ 

307 

308  wf  =  sqrt(  zo  +  (tau-tauw)*(tau-tauw)/zo  ); 

309  printf  ("w(“/„g)  =  “/og\n"  ,tau,  wf )  ; 

310 

311 

312  for(k=0;k<10;k++) 

313  for(l=0;l<10;l++) 

314  for(x=0;x<nx;x++) 

315  for(y=0;y<nx;y++){ 

316  u  =  (sqrt(2)/wf)*(x  -  nx/2  +0.5)  *  dx; 

317  V  =  (sqrt(2)/wf)*(y  -  nx/2  +0.5)  *  dx; 

318  if  (k==0)  Hx=1.0; 

319  if  (k==l)  Hx=2*u; 

320  if  (k==2)  Hx=4*pow(u,2)-2; 

321  if  (k==3)  Hx=8*pow(u,3)-12*u; 

322  if  (k==4)  Hx=16*pow(u,4)-48*pow(u,2)+12; 

323  if  (k==5)  Hx=32*pow(u,5)-160*pow(u,3)+120*u; 

324  if  (k==6)  Hx=64*pow(u,6)-480*pow(u,4)+720*pow(u,2)-120; 

325  if  (k==7)  Hx=128*pow(u,7)-1344*pow(u,5)+3360*pow(u,3)-1680*u; 
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326  if  (k==8)  Hx=256*pow(u,8)-3584*pow(u,6)+13440*pow(u,4)-13440* 

327  pow(u,2)+1680; 

328  if  (k==9)  Hx=512*pow(u,9)-9216*pow(u,7)+48384*pow(u,5)-80640* 

329  pow(u,3)+30240*u; 

330  if  (1==0)  Hy=1.0; 

331  if  (1==1)  Hy=2*v; 

332  if  (1==2)  Hy=4*pow(v,2)-2; 

333  if  (1==3)  Hy=8*pow(v,3)-12*v; 

334  if  (1==4)  Hy=16*pow(v,4)-48*pow(v,2)+12; 

335  if  (1==5)  Hy=32*pow(v,5)-160*pow(v,3)+120*v; 

336  if  (1==6)  Hy=64*pow(v,6)-480*pow(v,4)+720*pow(v,2)-120; 

337  if  (1==7)  Hy=128*pow(v,7)-1344*pow(v,5)+3360*pow(v,3)-1680*v; 

338  if  (1==8)  Hy=256*pow(v,8)-3584*pow(v,6)+13440*pow(v,4)-13440* 

339  pow(v,2)+1680; 

340  if  (1==9)  Hy=512*pow(v,9)-9216*pow(v,7)+48384*pow(v,5)-80640* 

341  pow(v,3)+30240*v; 

342 

343  rsq  =  ( (x-nx/2+0 . 5) * (x-nx/2+0 . 5)+ (y-nx/2+0 . 5) * (y-nx/2+0 . 5) ) *dx*dx; 

344  phi  =  (  rsq* (tau-tauw) / (wo*wo*wf *wf ) - (k+1+1) *atan( (tau-tauw) /zo)  ); 

345 

346  Abfinalr[k] [1] [x] [y]  =  Ao  *  wo  *  Hx  *  Hy  *  exp(-rsq/ (wf *wf ) )  * 

347  cos (phi) /wf *sqrt (1/ (pow (2 ,k) *pow(2 , 1) *f actorial (k) *f actorial (1) ) ) ; 

348 

349  Abf inali [k] [1] [x] [y]  =  Ao  *  wo  *  Hx  *  Hy  *  exp(-rsq/ (wf *wf ) )  * 

350  sin(phi)/wf *sqrt (1/ (pow(2 ,k) *pow(2 ,1) *f actorial (k) *f actorial (1) ) ) ; 

351 

352  Powerf inal [k] [1]  =  Powerf inal [k] [1]  + 

353  Abf  inalr  [k]  [1]  [x]  [y]  *Abf  inalr  [k]  [1]  [x]  [y]  + 

354  Abf  inali  [k]  [1]  [x]  [y]  *  Abf  inali  [k]  [1]  [x]  [y]  ;  } 

355 

356  printf  ("P_00(1)  =  "/ogXn"  ,  Powerf  inal  [0]  [0] )  ; 

357 

358  /*  Finding  total  power  in  At  at  tau=l  * (should  equal  the  tau=0  case)*/ 

359  f or (x=0 ; x<nx; X++) 

360  for(y=0;y<nx;y++) 

361  PI  =  PI  +  Atr[x]  [y]*Atr[x]  [y]  +  Ati  [x]  [y]  *Ati  [x]  [y]  ; 

362 

363  printf ("P_T(1)  =  “/.gXn",  PI); 

364 

365  /*  orthogonality  calculation  */ 

366  K  =  (2*dx*dx) / (pi*zo) ; 

367  for(k=0;k<10;k++) 
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368 

369 

370 

371 

372 

373 

374 

375 

376 


for(l=0;l<10;l++) 
for (x=0 ; x<nx ; x++) 
for(y=0;y<nx;y++)  { 
cfinalr[k]  [1]  =  cfinalr[k]  [1]  +  K  * 

(Atr  [x]  [y]  *Abf  inalr  [k]  [1]  [x]  [y]  +  Ati  [x]  [y]  *Abf  inali  [k]  [1]  [x]  [y] ) ; 


cf inali  [k] [1]  =  cf inali  [k]  [1]  +  K  * 

(Ati  [x]  [y]  *Abf  inalr  [k]  [1]  [x]  [y]  -  Atr  [x]  [y]  *Abf  inali  [k]  [1]  [x]  [y] ) ; } 


377  /*  calculating  cmn  squared  */ 

378  for(k=0;k<10;k++) 

379  for(l=0;l<10;l++)  { 

380  c2magfinal[k] [1]  =  (cf inalr [k] [1] *cf inalr [k] [1]  + 

381  cf inali [k] [1] *cf inali [k] [1] ) ; 

382  c2suinFinal  =  c2sumFinal  +  c2magf  inal  [k]  [1]  ; } 

383 

384  printf ("\nNormalized  initial  cmn  amplitudes  squared\n"); 

385  for(k=0;k<10;k++) 

386  for(l=0;l<10;l++)  { 

387  fS[k][l]  =  c2magfinal[k]  [l]/c2sumFinal; 

388  if  (fS[k][l]  >  0.01)  printf  ("Cmn2 ["/od]  ["/od]  =  y„g\n"  ,k,l,fS[k]  [1] )  ;} 

389  printf ("\n") ; 

390  printf ("\nClosing  f iles . . . \n") ; 

391  fclose (input) ; 

392  fclose(initial) ; 

393  fclose(f inal) ; 

394  fclose (ytau) ; 

395 


396  }  /*end  mainO  */ 


397 

398  int  factorial (int  a) 

399  { 

400  if  ((a==l)  II  (a==0)) 

401  return  1; 

402  else  { 

403  a  *=  factorial (a-1) ; 

404  return  a;  }  } 

405  /*  End  of  Cmn.c  */ 
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2.  CMN.IN  -  INPUT  FILE 

/*  Begin  input  file  */ 


1 

1  0 

2 

2 

1  2 

0 

3 

0  0 

2 

4 

0  4 

0 

5 

0  0 

4 

6 

0  2 

7 

7 

0  4 

9 

8 

0  3 

8 

9 

0  5 

3 

10 

0  0 

6 

11 

1.0 

10 

100 

12 

0.5 

1.1 

0  400E-6 

13 

y.co 

mO 

nO 

14 

7.C1 

ml 

nl 

15 

y.c2 

m2 

n2 

16 

y.c3 

m3 

n3 

17 

yoc4 

m4 

n4 

18 

y«c5 

m5 

n5 

19 

IcQ 

m6 

n6 

20 

y«c7 

m7 

n7 

21 

y.c8 

m8 

n8 

22 

y.c9 

m9 

n9 

23 

y.Ao 

W  nx 

24 

"/oZO 

tauw  dtau 

/*  End 

input  ; 

file  */ 

3.  COEFFS.M  -  MATLAB  SCRIPT 

1  close  all 

2  clear 

3  clc 

4  load  initial. out 

5  load  final. out 

6  load  Cmn.in 

7  load  ytau.out; 

8 
9 

10  ytau  =  ytau'; 

11  initial  =  initial'; 

12  final  =  final'; 
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13 

14 

15  W  =  Cmn(ll,2); 

16  nx  =  Cmn(ll,3); 

17  dtau  =  Cmn(12,3); 

18  tau=0 : dtau: 1-dtau; 

19  X  =  0 : 1 :nx-l ; 

20 

21  [npxjtimesteps]  =  size(ytau); 

22  dtau  =  1/timesteps; 

23  tau  =  0: dtau: 1-dtau; 

24  X  =  0:l:npx-l; 

25 

26 

27  [x_initial,y_initial]=meshgrid(0:W/nx:¥-¥/nx) ; 

28  [x_final,y_final]=meshgrid(0:¥/nx:¥-¥/nx) ; 

29 

30 

31  figure (' Color \ 'white 0 

32  setCgca, 'FontSize' ,25) 

33  pcolor (x_initial , y_initial , initial) 

34  colorbar 

35  axis  square 

36  shading  flat 

37  shading  interp 

38  title (' initial  waveform') 

39  xlabel('x') 

40  ylabel('y') 

41  colormap(hot) 

42 

43  figure ( 'Color' , 'white') 

44  setCgca, 'FontSize' ,25) 

45  surf (initial) 

46  title ('initial  surface  plot') 

47  xlabel('x') 

48  ylabel('y') 

49  zlabel('a_i*a_i') 

50  colormap(hot) 

51 

52  figure ( 'Color' , 'white') 

53  set(gca, 'FontSize' ,25) 

54  pcolor(x_final,y_final, final) 
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55  xlabeK^xO 

56  ylabeK^yO 

57  colorbar 

58  axis  square 

59  shading  flat 

60  shading  interp 

61  title (^final  waveformO 

62  colormap(hot) 

63 

64  figure (' Color \ 'white 0 

65  setCgca, 'FontSize' ,25) 

66  surf (final) 

67  xlabel('x') 

68  ylabel('y') 

69  zlabel('a_f*a_f*') 

70  title ('final  surface  plot') 

71  colormap(hot) 

72 

73  figure ( 'Color' , 'white') 

74  set(gca, 'FontSize' ,25) 

75  [Tau,X]=meshgrid(tau,x) ; 

76  pcolor (Tau,X,ytau) 

77  shading  flat 

78  shading  interp 

79  title ('Laser  amplitude  y  vs. 

80  xlabel('\tau') 

81  ylabel('x') 

82  colormap(hot) 


\tau' ) 
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APPENDIX  B.  USEFUL  RELATIONS 


Quick  reference  for  some  important  equations 


1.  FREE  ELECTRON  LASER  OVERVIEW 


Resonance  Condition 


Gain 


Extraction 


A  = 

G  = 

V  = 


P-Po 

Po 


_ Optical  power _ 

Initial  electron  beam  power 


Equation  (II.  1) 
Equation  (II. 2) 
Equation  (II. 3) 


2.  FREE  ELECTRON  LASER  THEORY 


Pendulum  Equation 


Power 


PEL  wave  equation 


c=  |a|  cos(C  +  (/>) 
P=|aP 


dr 


a(f,T) 


=<  —je  > 


Equation  (III. 32) 
Equation  (III. 37) 
Equation  (III. 64) 
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3.  OPTICAL  THEORY 


F.  Trial  Solution 

a  =  ttoS  V  ^09 ; 

Equation  (IV. 3) 

F.  Solution 

f"^  if  ’’y  — arctan 

a  =  tto—e 

^  W 

)  Equation  (IV.  20) 

HOM  Trial  Solution 

(  ^ 

a  = 

Equation  (IV.  24) 

HG  Mode  Solution 

am,n  = 

Equation  (IV. 45) 

HOM  Phase 

=  ^2^2  (m  +  n  +  l)  arctan  J  Equation(IV.46) 

Normalized  HG  mode 

/I  —  1  ^ 

^m,n  v'2™2"m!n! 

Equation  (IV. 60) 

Total  field  HG 

At  =  Y2m,n  Cm,nAm,n 

Equation  (IV. 61) 

Normalized  Goefficients 

Cm,n  =  ^ai^2  At  dxdy  Equation  (IV. 66) 

SIMULATIONS 

Total  Power  in  optical  field  Pt  =  Ylmn^m,n 

Equation  (V.13) 

Fractional  Power  Spectrum  S'[m][n]  = 

Equation  (V.15) 
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